530SW86006
EPA/530-SW-86-006
April 1986
SOILINER Model - Documentation
and User's Guide (Version 1)
Contract No. 68-01-6871
Assignment No. 48
In Cooperation with
Office of Solid Waste
U.S. Environmental Protection Agency
Washington, D.C. 20460
Hazardous Waste Engineering Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, Ohio 45268
-------
DISCLAIMER
This report was prepared by Russell A. Johnson, Eric S. Wood, and Richard
J. Wozmak of GCA Technology Division, Inc., Bedford, Massachusetts 01730,
under Contract No. 68-01-6871. The technical project monitors were Les Otte
of the Office of Solid Waste and Douglas Ammon of the Hazardous Waste
Engineering Research Laboratory. This is a draft report that is being
released by EPA for public comment on the accuracy and usefulness of the
information in it. The report has received extensive technical review but the
Agency's peer and administrative review process has not yet been completed.
Therefore it does not necessarily reflect the views or policies of the
Agency. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
11
-------
FOREWORD
The Enviromental Protection Agency was created because of increasing
public and governmental concern about the dangers of pollution to the health
and welfare of the American people. Noxious air, foul water, and spoiled land
are tragic testimony to the deterioration of our natural environment. The
complexity of the environment and the interplay of its components require a
concentrated and integrated attack on the problem.
Today's rapidly developing and changing technologies and industrial
products and practices frequently carry with them the increased generation of
solid and hazardous wastes. These materials, if improperly dealt with, can
threaten both public health and the environment. Abandoned waste sites and
accidental releases of toxic and hazardous substances to the environment also
have important environmental and public health implications. The Hazardous
Waste Engineering Research Laboratory assists in providing an authoritative
and defensible engineering basis for assessing and solving these problems.
Its products support the policies, programs, and regulations of the
Environmental Protection Agency; the permitting and other responsibilities of
State and local governments; and the needs of both large and small businesses
in handling their wastes responsibly and economically.
The Office of Solid Waste is responsible for issuing regulations and
guidelines on the proper treatment, storage, and disposal of hazardous wastes,
in order to protect human health and the environment from the potential harm
associated with improper management of these wastes. These regulations are
supplemented by guidance manuals and technical guidelines, in order to assist
the regulated community and facility designers to understand the scope of the
regulatory program. Publications like this one provide facility designers
with state-of-the-art information on design and performance evaluation
techniques.
This document describes technical procedures for determining adequate
thicknesses of compacted soil liners to prevent migration of hazardous
constituents through the liner during the active life and post-closure care
period. It includes a performance simulation model that is based on numerical
techniques.
Marcie Williams
Director, Office of Solid Waste
Thomas R. Hauser
Director, Hazardous Waste Engineering Research Laboratory
iii
-------
PREFACE
Subtitle C of the Resource Conservation and Recovery Act (RCRA) requires
the U.S. Environmental Protection Agency (EPA) to establish a Federal
hazardous waste management program. This program must ensure that hazardous
wastes are handled safely from generation until final dispostion. EPA issued
a series of hazardous waste regulations under Subtitle C of RCKA that are
published in 40 Code of Federal Regulations (CFR) 260 through 265 and 122
through 124.
Parts 264 and 265 of 40 CFR contain standards applicable to owners and
operators of all facilities that treat, store, or dispose of hazardous
wastes. Wastes are identified or listed as hazardous under 40 CFR Part 261.
Part 264 standards are implemented through permits issued by authorized States
or EPA according to 40 CFR Part 122 and Part 124 regulations. Land treatment,
storage, and disposal (LTSD) regulations in 40 CFR Part 264 issued on July 26,
1982, establish performance standards for hazardous waste landfills, surface
impoundments, land treatment units, and waste piles.
EPA is developing three types of documents for preparers and reviewers of
permit applications for hazardous waste LTSD facilities. These types include
RCRA Technical Guidance Documents, Permit Guidance Manuals, and Technical
Resource Documents (TRD's).
The RCRA Techical Guidance Documents present design and operating
specificatons or design evaluation techniques that generally comply with, or
demonstrate compliance with the Design and Operating Requirements and the
Closure and Post-Closure Requirements of Part 264.
The Permit Guidance Manuals are being developed to describe the permit
application information the Agency seeks, and to provide guidance to
applicants and permit writers in addressing information requirements. These
manuals will include a discussion of each set of specifications that must be
considered for inclusion in the permit.
The Technical Resource Documents present state-of-the-art summaries of
technologies and evaluation techniques determined by the Agency to constitute
good engineering designs, pratices, and procedures. They support the RCRA
Technical Guidance Documents and Permit Guidance Manuals in certain areas
(i.e., liner, leachate management, closure covers, and water balance) by
describing current technologies and methods for designing hazardous waste
facilties, or for evaluating the performance of a facility design. Although
emphasis is given to hazardous waste facilities, the information presented in
these TRD's may be used for designing and operating nonhazardous waste LTSD
IV
-------
facilities as well. Whereas the RCRA Technical Guidance Douments and Permit
Guidance Manuals are directly related to the regulations, the information in
these TRD's covers a broader perspective and should not be used to interpret
the requirements of the regulations.
This document is a draft being made available for public review and
comment. It has undergone review by recognized experts in the technical areas
covered, but Agency peer review processing has not yet been complete. Public
comment is desired on the accuracy and usefulness of the information presented
in this document. Comments received will be evaluated, and suggestions for
improvement will be incorporated, wherever feasible, before pulication of the
next edition. Communications should be addressed to Docket Clerk, Room
S-212(A), Office of Solid Waste (WH-562), U.S. Environmental Protection
Agency, 401 M Street, S.W., Washington, D.C., 20460. The document under
discussion should be identified by title and number; e.g., "SOILINEK Model -
Documentation and User's Guide (Version 1)" (EPA/530-SW-8b-006).
v
-------
ABSTRACT
This Technical Resource Document provides documentation and a User's
Guide for the SOILINER computer model. SOILINER is a finite-difference
approximation of the highly nonlinear, governing equation for one-dimensional,
unsaturated flow in the vertical dimension. SOILINER was designed to simulate
the dynamics of an infiltration event across a compacted soil liner system
beneath impounded liquid. Since the governing equation reflects liner
heterogeneity and the dependence of liner properties on the degree of
saturation, SOILINER is capable of accurately representing infiltration for a
variety of soil(clay)liner scenarios. Important features inherent to the
SOILINER model include the ability to simulate: (1) multilayered system, (2)
variable initial moisture content, and (3) changing conditions on the
boundaries of the compacted soil liner flow domain. Due to these features,
SOILINER serves as a comprehensive tool for the design of liner
configurations, specifically liner conductivity and thickness.
VI
-------
ACKNOWLEDEGMENTS
The authors would like to thank Dan Goode (now with the Nuclear
Regulatory Commission) and Dr. Michael Mills of GCA for their contributions to
the SOILINER Documentation and Users' Guide. SOILINER, a finite difference
model originally authored at GCA by Dan Goode, stood the test of time during
our extensive sensitivity analysis with very little modification. We were
continually surprised by the versatility of this model and the role it plays
as a tool in understanding the dynamics of an infiltration event. Also, the
original SOILINER document (Goode and Smith, 1984) was clear and precise, and
thus served as the basis for Part I of this report, where only minor changes
were made.
Dr. Mills provided technical guidance on many aspects of our most recent
effort to review and modify SOILINER. We acknowledge his invaluable
suggestions during development of new model features and his efforts to
periodically review the contents of this document. His accessibility and
eagerness to help is greatly appreciated and contributed immensely to the
production of the SOILINER Documentation and Users' Guide.
VII
-------
TABLE OF CONTENTS
Section
PART I - SOILINER DOCUMENTATION
1 INTRODUCTION 1-2
2 CONCEPTUAL MODEL OF FLOW THROUGH LINER 1-3
2.1 Description of Physical Problem 1-3
2.2 Mathematical Statements 1-3
2.2.1 Governing Equation of Vertical
Unsaturated Flow 1-3
2.3 Soil Moisture Characteristics 1-7
2.4 Unsaturated Hydraulic Conductivity 1-12
2.5 Characteristic Curves Available in SOILINER 1-16
3 NUMERICAL SIMULATION OF UNSATURATED FLOW 1-23
3.1 Finite Difference Method 1-23
3.1.1 FDM Spatial Difference Approximations ...... 1-24
3.1.2 FDM Temporal Derivative Approximations 1-24
3.2 Numerical Solution 1-27
3.2.1 Transient Finite Difference Equation 1-27
3.2.2 Steady State, Finite Difference Equation .... 1-30
3.3 Flux and Velocity Calculations 1-32
3.4 Method of Determining Breakthrough 1-33
4 VERIFICATION 1-34
4.1 No-Flow Steady State 1-34
4.2 Steady State Evaporation Flow Upward 1-35
4.3 Transient Infiltration 1-37
4.4 Particle Tracking 1-37
Vlll
-------
TABLE OF CONTENTS
Section Page
PART II - SOILINER USER'S GUIDE
1 INTRODUCTION II-2
1.1 History of SOILINER Computer Code II-2
1.2 Scope of the Manual II-2
1.3 Computer Requirements II-3
2 MODELING CAPABILITIES II-4
2.1 Modeling Assumptions II-4
2.2 Solution Strategies II-4
2.3 Definition of Clay Liner Domain II-5
2.4 Determining Breakthrough 11-10
2.5 Variable Initial and Boundary Conditions ....... 11-10
3 DATA REQUIREMENTS 11-15
3.1 Description of Input Data 11-15
3.2 Measurement of the Soil Moisture Characteristics . . . 11-17
3.3 Functional Relationship Formula 11-18
3,4 Characteristic Curves Available in SOILINER II-22
4 PROCEDURE FOR APPLYING SOILINER 11-26
4.1 Model Set-up 11-26
4.2 Designing a Finite Difference Grid 11-26
4.3 Steady State Solution Strategy 11-27
4.4 Transient Solution Strategy 11-29
4.4.1 Method of Temporal Sensitivity Analysis .... 11-34
4.4.2 Choosing a Time Step 11-35
4.4.3 Choosing a Temporal Weighting Parameter .... II-41
4.4.4 Setting a Tolerance for
Iteration Convergence 11-44
4.4.5 Choosing a Maximum Number of Iterations
Per Time Step 11-44
4.4.6 Sensitivity of the Particle Tracking Algorithm 11-46
REFERENCES 11-48
APPENDICES
A List of Symbols Used in Text A-l
B Gardner's Analytical Solution B-l
C Partial List of Variables Found in SOILINER C-l
D Example of Input Data D-l
E Example of Output Data E-l
F Subroutine Descriptions F-l
G SOILINER Source Code G-l
H SOILINER Package H-l
IX
-------
FIGURES
Number
PART I - SOILINER DOCUMENTATION
2-1 Flow domain for liner breakthrough 1-4
2-2 Schematic of initial capillary liquid pressure distribution . . 1-5
2-3 The effect of texture on soil moisture characteristics .... 1-8
2-4 A piecewise linear relation between moisture content
and the logarithm of suction 1-11
2-5a Sandy and clayey soil moisture characteristics computed
from a continuous functional relationship developed by
Milly and Eagleson (1980), moisture content 1-13
2-5b Sandy and clayey soil moisture characteristics computed
from a continuous functional relationship developed by
Milly and Eagleson (1980), log of specific moisture
capacity (cm~~l) vs. PF = log (-W, ^ in cm 1-13
2-6 Schematic of unsaturated hydraulic conductivity for a
sand and clay soil 1-15
2-7 The moisture characteristic curve using equations (2-27)
and(2-28) for the hyperbolic and parabolic sections
respectively (broken line segments are disregarded) 1-18
2-8a Characteristic moisture curves for soil types ranging
from sand to sandy clay loam 1-20
2-8b Characteristic moisture curves for soil types ranging
from silty clay loam to a compacted clay 1-21
3-1 Finite difference method spatial discretization using
mesh-centered grids 1-24
3-2 Finite difference method temporal discretization for
node i at time level n + 1 showing explicit (dash)
and implicit (solid) relationships 1-27
3-3 SOILINER solution procedure flow chart 1-31
4-1 Comparison of solutions for steady state vertical flow
upward from a water table 1-36
4-2 Comparison of a Phillips quasi-analytical solution and
SOILINER regular grid solution for infiltration into
Yolo light clay under ponding 1-38
4-3 Comparison of a Phillips quasi-analytical solution and
SOILINER graded grid solution for infiltration into Yolo
light clay under ponding 1-39
4-4 Changes in particle depth with time across a fully saturated
clay liner under steady state conditions .... 1-40
-------
FIGURES
Number Page
PART II - SOILINER USER'S GUIDE
2-1 Changes in pressure distribution for a clay liner
underlain by a high conductivity, native soil II-6
2-2 Changes in pressure distribution for a clay liner
underlain by a low conductivity, native soil II-7
2-3 Changes in pressure with time in the upper 300 cm
of a 482 cm clay-sand-clay liner configuration II-8
2-4 Changes in moisture content with time corresponding
to pressure curves shown in Figure 2-3 II-9
2-5 Changes in particle depth depicting the dynamic
characteristics of an infiltration event, specifically
decreasing wetting front velocities with time 11-11
2-6 Change in steady state pressure distribution due to
increased impoundment depth at year 1.5 11-13
2-7 Effect of increased impoundment depth on the time
to breakthrough 11-14
3-1 Illustration of the essential parts of a tensiometer 11-19
3-2 Illustration of the essential components of a portable
neutron soil-moisture meter 11-19
3-3 Qualitative illustration of a moisture characteristic
curve (broken line segments are disregarded) 11-20
3-4a Characteristic moisture curves for soil types ranging from
sand to clay loam 11-23
3-4b Characteristic moisture curves for soil types ranging from
silty clay loam to a compacted clay 11-24
4-1 Two vertical grids with variable subdomain sizes 11-28
4-2 Steady state algorithm leads to numerical oscillation
as shown at three selected nodes from a clay/sand
liner configuration (SRPARM = 1.0) 11-30
4-3 Impact of SRPARM on convergence during steady state
algorithm at node 85 11-31
4-4 SOILINER sensitivity to initial DT values (pressure
distributions at ENDTIM of transient simulations where
CHPARM = 75 and ALPHA = 1.0) 11-37
4-5 Divergence associated with large DT and CHPARM values .... H-39
4-6 Results of MAXIT warning due to excessive DTMAX, where the
magnitude of error between successive iterations (with
respect to ERRMAX) can be determined from ERR 11-42
4-7 Instability of particle movement due to inaccurate PSI
distributions associated with divergence at the time steps
shown above 11-47
XI
-------
TABLES
Number Page
PART I - SOILINER DOCUMENTATION
2-1 Representative Values of Hydraulic Parameters for a
Number of Soil Texture Classes 1-22
PART II - SOIUNER USERS' GUIDE
3-1 Representative Values of Hydraulic Parameters for the
Soil Texture Classes Coded into SOILINER II-25
4-1 Effect of SRPARM on the Number of Iterations to
Convergence (90 cm liner) ...... 11-32
4-2 Effect of SRPARM on the Number of Iterations to
Convergence (180 cm liner) 11-32
4-3 CPU Time as a Function of DT for Constant CHPARM = 100 .... 11-38
4-4 Sample of CPU Time Necessary for Simulations with
Small DT and CHPARM Values 11-40
4-5 Solution Characteristics as a Function of Temporal
Parameters and MAXIT 11-43
4-6 Values of DT and CHPARM Producing Register
Underflows/Overflows or Abends as a Function of Alpha .... H-45
Xll
-------
PART I
SOILINER DOCUMENTATION
1-1
-------
1. INTRODUCTION
On November 8, 1984 the Hazardous and Solid Waste Amendments of 1984
(HSWA) were signed into law. Among the provisions of this new law are minimum
technological requirements for hazardous waste landfills, surface
impoundments, and waste piles. These provisions are found in Sections 3004(o)
and 3015 of the Resource Conservation and Recovery Act (RCRA) and further
defined in the Part 264 regulation and associated guidance. Any new surface
impoundment must install two or more liners and a leachate collection system
between such liners. The lower liner must be designed, operated, and
constructed to prevent breakthrough of constituents over the period of
operation, including any post-closure care period.
Existing methods of evaluating flow through soil liners have been
reviewed, including the: (1) transit time equation, (2) Green-Ampt wetting
front model, and (3) transient, linearized, approximate model - see Appendix A
of the Technical Resource Document (TRD-EPA/530-SW-84-001). Simplifying
assumptions associated with these methods are fairly restrictive and do not
allow for an accurate simulation of the infiltration event.
SOILINER is a finite-difference approximation of the highly nonlinear,
governing equation of unsaturated, vertical flow. Since the governing
equation reflects liner heterogeneity and the dependence of liner properties
on the degree of saturation, SOILINER is capable of accurately representing
infiltration for a variety of compacted soil liner scenarios. Important
features inherent to the SOILINER model include the ability to simulate:
(1) multilayered systems, (2) variable initial moisture content, and
(3) changing conditions on the boundaries of the compacted soil liner flow
domain.
This report is comprised of two parts - documentation of the SOILINER
model and a User's Guide. Part I discusses the conceptual model of flow and
its mathematical representation (Section 2); the SOILINER finite-difference
formulation (Section 3); and model verification (Section 4). Part II presents
the modeling capabilities (Section 2), data requirements (Section 3), and
procedure for applying SOILINER (Section 4) with examples.
1-2
-------
2. CONCEPTUAL MODEL OF FLOW THROUGH LINER
2.1 Description of Physical Problem
The flow domain for liner breakthrough (shown in Figure 2-1) consists of
the following: a layer of liquid in the impoundment of depth h^ (L); a
compacted soil liner of thickness d (L); a layer of underlying site soil,
which may or may not be saturated, to a depth zw(L); and a constantly
saturated ground water layer of the same site soil. For this study, only
vertical processes are considered.
A schematic of the initial liquid pressure distribution (40 is shown in
Figure 2-2. Before installation of the liner, the soil moisture at the site
can be assumed to be in static equilibrium with the underlying water table and
saturated zone. Departures from this condition can occur if there is
significant evaporation from or recharge to the water table, and these
departures can be easily quantified. The soil liner is installed on top of
the site soil and is compacted. The liner is homogeneous and hence has an
initially constant moisture content and constant pressure over its entire
thickness.
After the impoundment is filled, the flow system is not in equilibrium,
and liquid will flow vertically down from the impoundment into the liner, and
eventually into the site soil and saturated ground water zone. Our goals are
to simulate this flow and to predict the liner thickness required to prevent
leachate from migrating through the liner into the subsoils during the active
life and post-closure care period of the unit.
2.2 Mathematical Statements
2.2.1 Governing Equation of Vertical Unsaturated Flow—
The governing equation for transient, unsaturated flow in the vertical
direction can be written:
= 0 (2-la)
0 L U If O'f. L rt Z J
and, for steady state:
|^ K(*) || + K(i|0~| = 0 (2-lb)
L -J
in which ^= 4>-z [L] is matric potential or capillary pressure head, where
[L] is piezometric head; 6 [-] is volumetric moisture content; K(40 =
Kr(40 Ks [LT ] is unsaturated hydraulic conductivity, where Kr(/3z)
and gravitational potential K(^). The flux terra is developed from the
generalization of Darcy's law for water flow in porous media:
1-3
-------
DATUM
J\
w
v
V
IMPOUNDED
LIQUID
/ SOIL
/ LINER
MATERIAL
/
/
'
NATURAL
FILL
^SATURATED,".
«;00« ZONE •» V*
oo0e>o«°p'
•a o O » *
Figure 2-1. Flow domain for liner breakthrough.
1-4
-------
V
IMPOUNDMENT
HYDROSTATIC
PRESSURE
A
SOIL
LINER
CONSTANT
INITIAL
PRESSURE
UNSATURATED
NATURAL
FILL
SATURATED
ZONE
HYDROSTATIC
PRESSURE
T
NEGATIVE (SUCTION) (POSITIVE
LIQUID PRESSURE
Figure 2-2. Schematic of initial capillary liquid pressure distribution.
1-5
-------
q = -KAY; r- = -KAY; ——— = -KAY; i T- -•- i i (2-2)
f) Z
in which q [LT~^] is flux. Some assumptions implicit in (2-1) are that the
fluid has constant density and does not freeze, that the medium does not
deform, that the air phase is always at a spatially constant atmospheric
pressure, and that the water flow is unaffected by temperature gradients or by
solute concentration gradients.
Equation (2-la), derived for unsaturated flow, is also applicable to
modeling temporarily or permanently saturated soil zones. In that case, vis
known as the pressure head, d6/d^ becomes the specific storativity, and K(ip)
is equal to Kg.
2.2.2 Boundary Conditions—
At the top of the liner, z=0, the pressure head is controlled by the
level of liquid in the impoundment (see Figure 2-2). In terms of piezometric
head, the matric potential at the liner top is fixed:
4* = 41 „ - z at z=0 (2-3)
z=0
Since the piezometric head is constant in space (hydrostatic) in the ponded
liquid,
and, since z=0, (2-3) becomes:
ty = ((> - o = h at z=0 (2-5)
z=0 i
This is the fixed pressure Dirichlet boundary condition applied at the top of
the liner.
By definition, the inatric potential at the water table is equal to zero;
thus at the column bottom, z = zw (where zw is a negative value - see
Figure 2-1), the Dirichlet boundary condition is:
= z or 4> = 0 atz = z (2-6)
w w
This water table boundary condition is assumed to be controlled by local
ground water flow and to be unaffected by the amount of liquid discharging
through the liner. That is, the water table elevation is assumed to be
constant.
The matric potential matching condition between the liner soil and the
site soil is [see Bear, 1979, p. 206]:
at z = -d (2-7)
s
where <[> and ^s, refer to the matric potential of the liner and site soil
respectively. This is simply a condition of pressure continuity. This
1-6
-------
condition is incorporated into the governing equation (2-la) and is implicitly
satisfied.
2.3 Soil Moisture Characteristics
Unsaturated flow is subject to negative pressures (also called suction)
and this potential constitutes a moving force (Hillel, 1971). The raatric
suction is due to capillary and adsorptive forces which tend to move water
from where the suction is low to where it is high.
In a saturated soil at equilibrium, the suction is zero at the free water
surface. Away from the free water surface, interstitial pore spaces are
capable of remaining saturated under slight suctions, an example of this being
the capillary fringe of a given soil. However, as suction is increased beyond
a critical value, the largest continuous pores begin to empty. This critical
suction is called the air entry pressure and is a function of soil type.
Typically, the air entry pressure is small in coarse-textured, well-aggregated
soils. As suction is further increased, more water is gradually taken out of
the soil as the larger pores drain first, followed by progressively smaller
pores. At some high suction value only the very narrow pores will retain
water. Increasing suction is thus associated with decreasing moisture
content. This function is represented graphically by a curve known as a
soil-moisture characteristic curve, which relates moisture content to pressure.
Figure 2-3 shows characteristic curves for a sandy and a clayey soil. It
has been observed experimentally that the soil moisture content matric
potential relationship is hysteretic; it has a different shape when soils are
wetting than when they are drying. If a soil sample was saturated at a
pressure greater than zero and the pressure was then lowered in a stepwise
fashion until it reached a level much less than zero, the moisture content at
each step would follow the drying curve. If water was added to the soil
sample, the pressure would follow the wetting curve. If the soil sample was
only partially wetted or dried, the pressure would follow internal lines which
are called scanning curves.
To evaluate the storage term of the governing equation (2-1), we must
define the soil specific moisture capacity C(^) [L~^]:
C(Ui) = $ (2-8)
which is the slope of the moisture characteristic curve in Figure 2-3, and is
calculated as the derivative of the moisture retention function:
0 = Q(4>). (2-9)
1-7
-------
on
o
Z
UJ
H-
O
O.
•CLAYEY SOIL
SANDY S01L-
- WATER CONTENT
Figure 2-3. The effect of texture on soil moisture characteristics
(after Hillel, 1971).
1-8
-------
Several mathematical expressions have been proposed for the moisture
retention function, or characteristic curve. Brooks and Corey [1964]
collected moisture retention data on many granular media and developed a power
law relationship for the drainage cycle:
n y; > iii,
-X b
r + (n " V f * - ^b (2'10)
b
in which 6r [-] is residual (nonreducible) moisture content, i/> 5 [L] is
bubbling or air entry pressure at which air first enters the draining column,
n [-] is porosity, and X is the pore size distribution coefficient. The
experimental data to develop this model included pressures of 0 ^ fy >_ -500
cm. The corresponding specific moisture capacity, which is the derivative of
is:
0 -A\ > ill,
D
(2-11)
'b
This function is discontinuous at saturation, 41 = ^. Clapp and Hornberger
(1978) designed a model using a power curve and a short parabolic section near
saturation to represent gradual air entry. This two part function was used to
relate soil moisture and raatric potential for a wetting soil. The power curve
representing the moisture characteristic for W < 0.92 is:
with soil wetness, W, equal to 6/6s where 6S is the saturated water
content. Both ips, the saturation matric potential, and the exponent b are
empirical and must be estimated. For 0.92 < W < 1, the moisture
characteristic can be described by the parabola:
where
d-w.)2 w.u-w.;
i 11
and
i mW.
i
is the matric potential at W = 0.92.
1-9
-------
King [1965] fit a hyperbolic function to moisture retention data for
several soil types with pressures as low as -100 cm:
n 6
cosh [ (
cosh [ (
+ 1) (2-14)
Lj O LJ
where ^15 = 1.5 x 10^ cm (15 bars) and 0^5 is the moisture content when
fy = ^15• Again, the derivative of 0 is discontinuous at saturation ( "
^Ij) but only requires three measured parameters: n, i^, and 815. This
model performed well for several soil types, including a clay.
McQueen and Miller [1974] proposed a linear relationship between pF =
log (-^), with ty in cm, and moisture content. The three straight segments
were a capillary segment from saturation to pF 2.5, an adsorbed segment from
pF 2.5 to 5.0 and a tightly adsorbed segment from pF 5.0 to 7.0. As pointed
out by McQueen and Miller [1974], their fitting procedure allows convenient
approximation of the entire moisture characteristic curve from few data points.
Milly and Eagleson [1980] developed a related continuous function,
considering only two unsaturated segments:
8(pF) = ^ In [exp |M (aT -s, pF)> + exp { M (a - s pF)>]
M ( L L J ( Z Z J
(2-15)
- rrr ln [exp + exp
N ( / / J
in which a , a , s , s and 0 are defined by Figure 2-4, and M and M1
control the curvature of the joining segments. On a linear portion of the cap-
illary range between pFm£n and pFg, this function can be approximated as:
= a2 - s2 log (-*) (2-16)
and the specific moisture capacity is approximately
1-10
-------
0, -^_
saturation capillary
range range
adsorption
range
Figure 2-4.
A piecewise linear relation between moisture content and the
logarithm of suction (after Milly and Eagleson, 1980).
1-11
-------
Since (2-15) has continuous slope, the exact form of (2-17), for all \p , is
also continuous and can be written:
.V1
(2-18)
_
M1
where:
M a,
P! - e
M a.
M'a,
P3 =
M1
P4 " e
u
P3 (-40
q = - M s log (e)
q2 = - M
log (e)
(2-19)
= -M's2 log (e)
Equations (2-15) and (2-18) are plotted for sandy and clayey soils in
Figure 2-5.
2.4 Unsaturated Hydraulic Conductivity
Darcy's Law for unsaturated flow is the same for saturated flow with the
provision that K is a function of suction (negative pressure). The equation
states that the flux rate is proportional to the total piezometric head
gradient times the hydraulic conductivity (Bear, 1979):
a*
(2-20)
For saturated flow, hydraulic conductivity is a function both of the soil
type and of the fluid. In this analysis, hydraulic conductivity will refer to
the value applicable when the fluid is the proposed impoundment leachate.
Furthermore, it will be assumed that the properties (density and viscosity) of
the leachate do not differ significantly from those of the native water.
For a given soil, the unsaturated hydraulic conductivity is dependent on
the moisture content or matric potential (suction) in the soil. As a function
of moisture content, 6, the hydraulic conductivity function shows little
hysteresis [Bear, 1979; Maulem, 1976] and since, as stated previously, S(ip)
hysteresis for this infiltration problem is ignored, the unsaturated hydraulic
1-12
-------
a.
15 .20 .25 .30 .35 .40
MOISTURE CONTENT
b.
1.5 2.0 2.5 3.0
Figure 2-5.
Sandy and clayey soil moisture characteristics computed from
a continuous functional relationship developed by Milly and
Eagelson (1980), a. moisture content and b. log of specific
moisture capacity (cm ) vs. PF = log(-i^), t^in cm.
1-13
-------
conductivity as a function of matric potential can also be determined
uniquely. Figure 2-6 shows a schematic for typical relationships between
hydraulic conductivity and matric potential for a sandy and a clayey soil.
As can be seen in Figure 2-6, when the pores are saturated (1^5 < ^ ^0)
the hydraulic conductivity is maximal. However, when the soils begin to drain
(ty < ^b) , some of the pores become air-filled and the conductive portion of
the soil's cross-sectional area decreases correspondingly. In addition, as
suction increases, the first pores to empty are the largest ones, which are
the most conductive. Thus, the water is redirected to the smaller, less
conductive pores.
It is well understood that under saturated conditions, a sandy soil
conducts water more rapidly than a clayey soil. However, under unsaturated
conditions, the opposite may occur. Since a sandy soil contains relatively
larger pores, these pores quickly empty and become non-conductive as suction
increases, thus, steeply decreasing the initially high conductivity value
associated with sandy soils. On the other hand, in a soil with smaller pores
such as clays, many of the pores remain full and conductive even at
appreciable suction. Thus, the conductivity does not decrease as fast in a
soil with small pores, and may actually become greater than that of a soil
with large pores subject to the same suction. In the field where unsaturated
conditions exist, it often happens that water flow is greater and persists
longer in clayey soils than in sandy soils.
Numerous functional relationships have been proposed for the unsaturated
hydraulic conductivity. These models include
Gardner [1958]: K(ip) = a / (b + (-Wm) (2-21)
where a, b, and m are constants;
Gardner [1958]: K(40 = Kg exp (-aty) (2-22)
where Ks is the saturated hydraulic conductivity;
Brooks and Corey [1964]: K(^) = K
a
(2-23)
for i|J <_ >Jj|j, where ipj, is the air entry pressure and X is an index of
pore-size distribution, and
Maulem [1978]: K(40 = K (S )0'015 w + 3'° (2-24)
S 6
1-14
-------
Ksat
o»
o
Figure 2-6. Schematic of unsaturated hydraulic conductivity for a
sand and clay soil.
1-15
-------
where Se is effective saturation:
S =
e n -
and
w
-
f Y i/i d9 (2-26)
ij; = +oo W
in which Yw is the specific weight of water [F/L3]. These formulas (2-24
- 2-26) require use of cgs units. Equation (2-26) "represents the amount of
work required to drain a unit volume of a saturated soil" [Maulera, 1978].
Thus, equation (2-24) is dependent on the moisture characteristic curve. This
model was shown to improve the prediction of unsaturated hydraulic
conductivity for fine-grained soil [Maulem, 1978] relative to previous
techniques.
There are several other integral techniques which derive the hydraulic
conductivity function from the moisture characteristic curve [see Maulem 1976;
Elzeftawy and Cartwright 1981; Jackson et al. , 1965]. Such techniques are
especially attractive when used in conjunction with numerical models because
of the numerical integration required. The accuracy of these methods, of
course, depends on the accuracy of the soil moisture characteristic curve, and
on the applicability of the physical assumptions made in a particular
technique.
2.5 Characteristic Curves Available in SOILINER
During the numerical solution of the governing equation for vertical
unsaturated flow, the soil moisture content, 6, and the hydraulic
conductivity, K, must be determined. The characteristic curves used by
SOILINER to relate soil moisture content with matric potential, fyt are modeled
using a power curve function combined with a parabolic function near
saturation to represent gradual air entry. This method is described in detail
by R.B. Clapp and G.M. Hornberger (1978). A conductivity function presented
by Campbell (1974) is used in SOILINER to relate soil moisture content with
unsaturated hydraulic conductivity, K. This function is estimated from the
soil moisture power curve stated above and the soil's saturated hydraulic
conductivity.
The power curve function representing the moisture characteristic is
* = *s W"b (2-27)
with soil wetness, W, equal to 9/6s, where 6g is the saturated water
content, or porosity. Both the matric potential at saturation, 4>8, and the
exponent b, which is the slope of the line segment for W < W^, are empirical
and must be estimated using a power curve regression analysis performed on the
log transforms of the parameters ^ and W. It has been shown that b is a
function of soil texture, and increases as soil texture moves from coarse to
1-16
-------
fine grained. The use of this function implies a sharp discontinuity in
matric potential near saturation. This is illustrated in Figure 2-7 as the
dotted line segment for W > Wj_. Although it is thought that some coarse
grained sands may have a small matric potential at W = 1, most soils,
particularly medium to finer grained soils show a gradual air entry region
near saturation. A modification to equation 2-27 is made to account for this
gradual air entry. Along a typical characteristic curve, there exists a point
where d^/dW changes from an increasing function to a decreasing function as W
increases. This inflection point has the coordinates (W^, 4^) as shown in
Figure 2-7, and the interval W^ < W < 1 is described by the parabola:
* = -m(W-n)(W-l) (2-28)
The parameters m and n are calculated such that:
1. equation 2-28 passes through points (W^, 4^) and (1,U) and,
2. d4>/dW of equations 2-27 and 2-28 are equal at the inflection point
The expression for the parameters m and n are
\\> . 41 .b
i i
U-wJ2
(W. (1-W.)
m =
(2-29)
i i
n = 2W. - -- -1
i mW.
i
Campbell (1974) derived a simple formula for estimating the unsaturated
conductivity given the matric potential at a specific time. This equation is
k = W2b+3 (2-30)
where k = K/KS. This formula has proven to be reasonably accurate over a wide
range of b values (wide range of soil textures) and for W values near
saturation.
In order to utilize equations 2-27, 2-28, and 2-30 for determining 6 and
Ks, for a given soil and matric potential, values for b, 4>8, 08, Ks, m
and n must be known. Table 2-1 provides representative average values of
these parameters for 11 soil textures. These values were obtained by applying
the power curve to desorption data from 1446 soil types. The soil samples
were collected from 34 localities throughout the United States. Each soil
type was assigned to a textural class based on its clay content. Values for m
and n corresponding to each textural class were estimated by setting W^ to
0.92. According to Clapp and Hornberger (1978), there was nothing in the
soils data indicating where this inflection point may occur. Rogowski (1971)
stated that this point usually occurs in the interval 0.8
-------
Z
Ul
H
O
a.
o
oc
s
Wi 1.0
SOIL WETNESS W
Figure 2-7.
The moisture characteristic curve using equations 2-27 and
2-28 for the hyberbolic and parabolic sections, respectively
(the broken line segments are disregarded).
1-18
-------
positive. For m to be positive, V^ must be greater than b/(b+l).
Figures 2-8a and 2-8b show the results of the two part function using the
parameters in Table 2-1 and the appropriate values for m and n.
The use of these equations in SOILINER makes it possible to determine,
for each value of matric potential, the hydraulic conductivity and moisture
content at each node. The advantages of selecting this method as opposed to
other methods described in Sections 2.3 and 2.4 are as follows:
1. With this approach only four parameters are needed (b, ^s, 6S>
Ks) to give a description of the hydraulic properties of a soil.
In addition, this method accounts for gradual air entry near
saturation.
2. The parameter, b, can be estimated from matric potential/moisture
content data or taken from Table 2-1 and used in equation (2-30), so
that unsaturated hydraulic conductivity need not be measured
directly.
3. The equations are straightforward and can be used to derive other
functional relationships if accurate soil moisture data is available.
4. Since this method categorizes soils into textural groups and uses an
average characteristic curve for each group, the user can simulate
soils in which no moisture suction data is available. This is done
by simply choosing the appropriate textural class for that soil.
The disadvantages of using this method to relate ty, Q , and K in the
SOILINER model include the following:
1. The values presented in Table 2-1 were derived from soil moisture
data obtained during the desorption process. Since the model
usually simulates a wetting front moving through a liner,
inaccuracies (due to the effects of hysteresis) may occur.
2. Errors may occur in using average values for each textural class as
opposed to values pertaining specifically to an individual soil and
its associated characteristic curve.
The equations have been useful to Clapp and Hornberger in several
different soil moisture models, one of which estimated the wetting front
suction required by the Green-Ampt equation. However, the average values
presented in Table 2-1 have not been verified and should be used with this
limitation in mind. In addition, these average parameters should be used only
in the absence of soil moisture data specific to a site. If soil moisture
data is available or can be obtained, it is recommended that the above
described functions be fitted to these data to obtain the most accurate
results. Sections 3.2 and 3.3 of the SOILINER User's Guide provides methods
for obtaining soil moisture data and describes how these functions are applied
to soil moisture data.
1-19
-------
I
Ni
o
C
0)
•*-•
C
o
o
ID
o
2
0)
1>
0.9 - SOIL TYPE 1
0
PF=(LOG(-PSI))
Figure 2-8a. Characteristic moisture curves for soil types ranging
from sand to sandy clay loam.
-------
I
S3
c
0)
c
o
o
0)
3
-*-»
OT
'o
2
0)
^
0)
0.1 -
0
PF=(LOG(-PSI))
Figure 2-8b.
Characteristic moisture curves for soil types ranging from silty clay loam to
a compacted clay.
-------
TABLE 2-1. REPRESENTATIVE VALUES OF HYDRAULIC PARAMETERS
FOK A NUMBER OF SOT1, TEXTURE CLASSES
Soil Mean clay
texture fraction
Sand
Loamy sand
Sandy Loam
Silt loam
Loam
Sandy clay loam
Silty clay loam
Clay loam
Sandy clay
Silty clay
Clay
0.03
0.06
0.09
0.14
0.19
0.28
0.34
0.34
0.43
0.49
0.63
b
4.05
4.38
4.90
5.30
5.39
7.12
7.75
8.52
10.4
10.4
11.4
*8
(cm)
12.1
9.0
21.8
78.6
47.8
29.9
35.6
63.0
15.3
49.0
40.5
e
g
(cm^/cm^)
0.395
0.410
0.435
0.485
0.451
0.420
0.477
0.476
0.426
0.492
0.482
K *
s
(cm/min)
1.056
0.938
0.208
0.0432
0.0417
0.0378
0.0102
0.0147
0.0130
0.0062
0.0077
* From Li et al. (1976)
1-22
-------
3. NUMERICAL SIMULATION OF UNSATURATED FLOW
Due to the nonlinearities of the unsaturated flow equation, exact
analytical solutions have not been obtained except for a few simple cases.
Numerical techniques are available to solve the unsaturated flow equations.
These techniques, primarily finite difference and finite element methods,
provide solutions of the governing equations that take into account the
inherent nonlinearities of the system and the heterogeneity of soil types and
variable initial conditions. They all rely on the basic concept of solving
for the moisture state at a finite number of points in space and time.
SOILINER utilizes the Finite Difference Method (FDM) to solve the
nonlinear, unsaturated flow equation. The FDM has been applied to vertical
unsaturated flow by many authors. Freeze [1969] investigated natural ground
water recharge and discharge mechanisms, firutsaert [1971] used a
two-dimensional vertical model in a study of soil moisture flow beneath drains
and irrigation ditches. Cooley [1971] investigated flow to a pumping well
using an axisymmetric two-dimensional vertical model. For one-dimensional
vertical flow, the FDM has been verified by, among others, Green et al.
[1970], Ragab et al. [1982], Giesel et al. [1973], and Elzeftawy and Dempsey
[1976]. Kunze and Nielson [1982], Haverkamp et al. [1977], and Reeder et al.
[1980] show excellent agreement between finite difference solutions and
Philip's [1958, 1969] quasi-analytical results. Unlike previous computer
models mentioned above, SOIUNER is specifically designed to evaluate soil
liner performance. Although SOILINER can be used in other applications it is
best utilized to simulate infiltrating leachate through soil profiles beneath
impounded liquid.
3.1 Finite Difference Method
The FDM consists of first breaking the solution domain of a differential
equation into subdomains. This process is called discretization. For each
subdomain, continuous differential terras in the governing equation are
replaced by approximate expressions based on the values of the state variable
in the subdomains. Typically, one or more equations are solved for each
subdomain. This technique is relatively simple to understand in practice, and
is applied to both spatial and temporal differential terras. Freeze and Cherry
[1979] describe the application of finite difference techniques to vertical
unsaturated flow.
In application to the vertical unsaturated flow equation, the vertical
column is divided into a row of short vertical segments. This row of segments
is called a grid. SOILINER utilizes mesh-centered grids, where the values of
matric potential (fy) are evaluated at nodes, which are located at the ends of
each segment. The vertical segments, or elements, may exhibit either variable
(Figure 3-la) or constant length (see Figure 3-lb). Likewise, solutions in
time are obtained by breaking time into discrete steps. Solutions at new
times are obtained (at the nodes) using the previous solution(s). In general,
the accuracy of this method improves as the grid spacing (Az) and the time
step (At) decrease in size.
1-23
-------
o. VARIABLE SPACING
b. CONSTANT SPACING
TOP
OF DOMAIN'
NODE 1
NODE 2
SOIL 1
SOIL
BOUNDARY
SOIL 2
NODE
NUMEL
BOTTOM
OF
DOMAIN
NODE
NUMEL + 1
SUBDOMAIN 1
SUBDOMAIN 2
SUBDOMAIN
NUMEL
NUMEL = number of subdomoins
or elements
Figure 3-1. Finite difference method spacial discretization
using mesh-centered grids.
1-24
-------
3.1.1 FDM Spatial Difference Approximations—
The spatial domain of the vertical unsaturated flow equation (2-la) is
the soil column, from the top of the impoundment liner down to the water
table. This domain includes the liner and the underlying site soil, and can
be discretized as shown in Figure 3-1 where NUMEL + 1 represents the bottom
node number. Values of the state variable (matric potential, ^), and soil
properties (hydraulic conductivity, K(i[>), and moisture capacity, C(<|0) are
approximated by values at a node i: ^£, K£, and G£. Note that K£ and
G£ are functions of ^£. These values are used to approximate the
derivatives in (2-la). The first term in (2-la) representing pressure driven
flux, can be replaced by:
Az I Az
in which z = zi+l~zi is tne constant node spacing and
T, f -ir ___ IT \ *• I £-
,1/2 »-2)
are the geometric mean hydraulic conductivities between nodes. The
gravitational flux term can be expressed as:
3K(i|) = i+l " i-l (3-3)
9z- 2Az
Equations (3-1) and (3-3) are centered difference approximations [see Finder
and Gray, 1977]. These expressions assume that the node spacing, Az, is
constant, although that is not a general requirement (Section 4 presents a
variable spacing, centered-difference algorithm analogous to 3-1). With
constant node spacing (3-3) is second order accurate. This means that if
conductivity K(IJJ) is a linear function of z and z^ only, then (3-3) is exact.
3.1.2 FDM Temporal Derivative Approximations—
The temporal domain of the vertical unsaturated flow equation (2-la) is
time, t > 0, after infiltration into the liner begins. 'Time is broken down
into time level subdomains in which superscript n represents time level. The
nodal values of properties and state variables, which are actually continuous
in time, are approximated by values at discrete time levels: 4"1 £, Kni,
and Cn £. To maintain computational efficiency, time level subdomains
(At) are automatically varied within SOILINER. First, an initial At and time
step change parameter are set, where At is multiplied by the change
parameter. Alone, the product of these two parameters would result in an
ever-increasing At. However, for each successive time step, At is also
modified by the maximum change in matric potential (") at any given node in
the flow domain between the previous and current time levels.
1-25
-------
Another feature incorporated into the SOILINER temporal derivative
approximation is that the storage term in (2-la) can be written:
C
3t n+a
n At
- C n+a (3-4)
in which Cn+a = a cn+1 + (1-a) Cn; 0 <_ a <_ 1 is a temporal weighting
parameter, and t = t11"1"^ - tn is the time step size. Subscripts 'have been
dropped for convenience. If a = 0, (3-4) becomes
At
which is forward differencing. When solving for ^n+^, all other terms of
(3-5) are known from the last time step and thus 4
-------
TIME LEVEL n
TIME LEVEL n +1
TIME LEVEL n+2
NODE i - 1
NODE i
NODE i +
KEY:
> IMPLICIT
•> EXPLICIT
Figure 3-2. Finite difference method temporal discretization for
node i at time level n + 1 showing explicit (dash) and
implicit (solid) relationships.
1-27
-------
in which C(^) = 36/8iJj [L"1] is specific moisture capacity, where 6 [-] is
volumetric moisture content; \\i [L] is raatric potential or capillary pressure
head; K(^) [LT"1-] is vertical unsaturated hydraulic conductivity; z [L] is
the vertical coordinate, negative downward; and t [T] is time. For notational
purposes, the second term in (3-7) which is flux divergence, is denoted by:
and by (3-7)
V(i|0 = C(<|0 f£ (3-9)
a<-
thus V [T~l] can be thought of as the rate of storage change at any point.
Standard centered finite difference expressions can be used to evaluate V
in (3-8). The spatial domain is discretized into a mesh centered grid (see
Figure 3-1) with node points on the boundaries between soil elements.
Substituting the FDM representations for the two flux terms (see (3-1) and
(3-3)) into (3-8), we have:
.
Az.+
(3-10)
Az.
i
where subscript "i" designates node number and
A2i+l/2 = zi+l - zi
Azi-l/2 = zi - zi-l (
Az. - 1/2 (z.+1 - Vl>
are the lengths of the two elements on each side of node i and the length
associated with node i, respectively. The element conductivities are
evaluated as geometric averages:
Ki+l/2 =
(3-12)
Ml/2
A weighted temporal difference expression is obtained by substituting
(3-9) into (3-4):
1-28
-------
in which 0 •" « ' 1 is the weighting parameter; and
_n+ot
i- (1-cO Cn (3-14)
These expressions apply at each node and subscripts have been dropped for
convenience.
Equation (3-13) is a FDM expression of the governing equation (3-7) at a
node. Since both Cn+1 and Kn+1 depend on the solution, i^n+I, equation
(3-13) is nonlinear in . Thus, after the initial prediction of ^n+l
(equation 3-13), an iterative procedure is used to resolve ^n+^. The
solution at a new iteration can be expressed as the solution from the last
iteration plus a correction computed at the new iteration
(3-15)
in which all terms are at node i and time step n+1 and superscript k is
iteration level. The left hand side (LHS) of (3-13) becomes:
At
At
(3-16)
The second term on the right hand side (RHS) of (3-13) is known from the last
time step. The first RHS term is evaluated by substituting (3-15) into (3-10)
v . v -
"i+l/2\
i
AJ+A***1-*^
Li-l/2\ A.._1/2
Azi+l/2
* ,k+1\
- A*. A
1 ,
Az.
(3-17)
1-1/2
or
v.
Azi
(3-18)
in which V(i|> ) is evaluated from the last iteration, and, again, superscript
"k" or "k+1" represents iterative values at time step n+1.
Tne final form of the FDM governing equation is obtained by substituting
(3-16) and (3-18) into (3-13) and grouping
1-29
k+1
terms on the LHS:
-------
(3-19)
All terms on the RHS of (3-19) are known from the last time step, "n", or the
last iteration, "k" at the new time step "n+1".
Equation (3-19) is repeated for each node of unknown matric potential
(i.e. excluding boundary nodes of fixed ty ). Thus, application of FDM to the
governing flow equation results in a system of equations to be solved each
time step for the unknown nodal values of matric potential. The resulting
matrix equation is solved directly forA^k+^ using the Thomas algorithm for
tridiagonal matrices (see Finder and Gray, 1977). Equation (3-19) is solved
iteratively for A^ at each time step, updating ty as shown in (3-15),
ultimately resolving the equations until the soil properties for that time
step approach constant values, andA^->0. Convergence at a given time step is
determined by a set criterion for the maximum Ai|> at any node between
successive iterations.
SOILINER will continue to iterate until the set criterion is achieved or
the maximum number of iterations specified per time step is exceeded, and a
forced exit from the iterative procedure occurs. In some cases convergence
can be achieved by simply increasing the number of maximum iterations, however
forced exits generally indicate divergence (i.e., A^->-a>). Conversely,
convergence may not occur if the set criterion is too stringent. For example,
the largest calculated Aij; at any node between successive iterations may
oscillate around a value of 0.01 cm but not converge if the error criterion
were set at 0.001 cm. Thus, convergence is a relative term associated with
the desired level of accuracy and the effort required to optimize a
combination of parameters for the solution strategy. Part II of this document
provides guidelines for the development of an appropriate set of parameters
for both transient and steady state solutions.
Figure 3-3 provides a flow chart for the overall procedure of a transient
solution strategy. When convergence occurs at a given time step, flux and
velocity calculations are made (Section 3.3), a particle is advected through
the liner as a means of determining breakthrough (Section 3.4), and a new time
step is determined. The procedure continues until steady state is achieved,
at which time the particle is tracked until the point of breakthrough.
3.2.2 Steady State, Finite Difference Equation—
A finite difference approximation of (2-lb) for steady state conditions
is developed identically to (3-19). However, the storage term C(^) and a are
set equal to 0.0 and 1.0 respectively in (3-19), resulting in:
Az.
1-1/2
k+ 1
. k+ 1
Azi-l/2
(3-20)
1-30
-------
NO
NliW TIME STEP n+1
COMPUTE V0|in)
SAVE if n
t
CALCULATE NEW TIME STEP
At
t
EXTRAPOLATE TO NEW
^n+1 USING V(i|/n)
COMPUTE FINAL
SOIL PROPERTIES
t
TRACK PARTICLE
T
CHECK MASS BALANCE
t
OUTPUT SOLUTION
IF REQUESTED
i
/^TEADY STATE^^
\ACHIEVED? /
YES
1
COMPUTE SOIL PROPERTIES
*" FOR*
t
COMPUTE V(i)ik)
t
BUILD STIFFNESS MATRIX
LHS OF GOVERNING EQUATION
Y
I
BUILD FORCING VECTOR
RHS OF GOVERNING EQUATION
*
SOLVE GOVERNING EQUATION
USING THOMAS ALGORITHM
1
\
YES ./HAS ^s. NO
>/ SOLUTION \v
\CONVERGED? >/
1—NO CBREAKTHROUGH?
Figure 3-3. SOILINER solution procedure flow chart,
1-31
-------
Since tli, in some cases it is still necessary
to use a less stringent convergence criteria than that which works for the
transient solution.
3.3 Flux and Velocity Calculations
Once the pressure distribution has been determined at any given time
step, or at steady state, flux and velocities are calculated across each
element in the flow domain. Element fluxes and velocities were chosen for two
reasons: (1) soil properties are assigned by element, and (2) a potential
gradient across each element can be determined by the calculated <|» values at
the two nodes which define a given element.
The flux term developed from a generalization of Darcy's law (2-2) is
approximated numerically as follows:
i+1
Az
+ 1.0
(3-22)
where qe is the element flux [LT l], Ke(^) is the interblock,
geometric-mean conductivity [LT'1], ^i+L and * £ are the potentials at
the bottom and top nodes of a given element respectively [L], and Az is the
element thickness [L]. Note that a negative value for qe indicates downward
flux, whereas a positive value indicates movement up, into the liner.
Equation (3-22) can be written as:
(3-23)
Az
where the first term on the RHS of (3-23) represents flux due to the matric
potential gradient and the second term represents flux due to gravitational
potential. Upon reaching steady state, all element fluxes should be equal
thus satisfying (2-lb).
1-32
-------
A conservative approach, with respect to liner design life, is provided
in SOILINER for element velocity calculations:
ve = qe/6e (3-24)
where ve is the interstitial pore velocity for a given element, qe is the
element flux, and 0g is the average moisture content for the nodes defining
the element of concern. Calculations .involving 6e, as opposed to porosity
(n), result in higher element velocities since the cross-sectional area of a
partially saturated pore capable of transmitting fluid is less than the actual
pore size (i.e., 0 distribution (and thus &
distribution) is achieved sooner than a particle can be advected across the
liner.
For the reasons stated above, a particle-tracking algorithm has been
included in the SOILINER model. This algorithm bases particle movement on
advection only and does not take into account the effects of dispersion.
However, element velocities used to advect the particle are a function of
space and time (i.e. ve = f(z,t))and thus reflect the dynamics of decreasing
velocities over time near the wetting front.
The particle is initially positioned at the liner surface and a velocity
for the first element is determined at to =0. A new ty distribution is then
calculated at t^ = t + At (as discussed in Section 3.2) from which a
second element velocity is calculated. Since the time-integration scheme
employs Simpson's rule, another time step is completed before particle
movement. At the time of particle movement, SOILINER then numerically
integrates the velocity function to determine the distance traveled over the
time period integrated. Finally, the total distance traveled over time is
updated and the particle is "located" with respect to the grid geometry such
that new element data is available for velocity calculations. Breakthrough
occurs when the particle passes a specified node point, generally set at the
liner/underlying-soil interface.
1-33
-------
4. VERIFICATION
In order to verify the numerical technique used by SOILINER, the model's
results were compared to analytical solutions for three physically realistic
problems. The first two simulations are for steady state unsaturated flow and
the third is for infiltration under ponding into an unsaturated clay.
Two grids were used for SOILINER's simulations. For each grid, node 1 is
at the column top, z = 0, and the last node is at the column bottom,
z = -50 cm. The first grid divided this 50 cm column into 50 1 cm elements
(51 node points). The second grid divided the soil column as follows: the
first 5 cm was divided into 10 elements; the next 25 cm into 25 elements; and
the final 20 cm into 10 elements for a total of 45 elements and 46 nodes.
4.1 No-Flow Steady State
The governing equation for steady state unsaturated flow is:
(4"u
which states that the flux is constant in space. The pressure boundary
conditions to (4-1) determine the direction and rate of flow.
At steady state no-flow, the piezometric gradient is zero:
i* = W*«) = o (4-2)
3z 8z
Thus, the capillary pressure head gradient is equal to -1:
li.-|«.-l (4-3)
3z 3z
The boundary conditions for this situation are
^ = 0 z = z (4-4)
w
at the water table at the bottom of the soil column and
z - z z=z (4-5)
w t t
at zt, the soil column top. Between the column bottom and top, capillary
pressure head is inversely proportional to elevation:
- z z < z < zt (4-6)
w w — — t
Equation (4-6) is an analytical solution to (4-1) with boundary conditions
(4-4) and (4-5).
1-34
-------
For the SOILINER simulations, the top node (zt = Ocm) capillary
pressure head is ^ = -50 cm and the bottom node (zw =-50cm) capillary
pressure is ^ = 0 cm. Simulations with both constant node spacing and
variable node spacing produce the exact analytical solution to five
significant digits.
4.2 Steady State Evaporation Flow Upward
Gardner [1958] developed analytical solutions for steady state
unsaturated flow when the unsaturated hydraulic conductivity follows the
following form:
(4-7)
where a and b are constants, and n is 1, 1/2, 2, 3, or 4. For n = 4, the
analytical solution to (4-1) is given implicitly by:
In
(\ii — p iK/2 + P \
__ _ _ I
_ / p^v/5\
V"wj.
(4-8)
in which
r = q/a
p4 - B/r
6 = rb + 1
where q is the discharge rate, constant in space and time, and W is a constant
of integration. For a water table at an elevation of zw = -50 cm and flow
upward, W = -50 cm.
Using (4-7) for the functional soil hydraulic conductivity, SOILINER was
run with boundary conditions of
and
= 0 at zw = -50 cm
= -200 cm at z = 0
(4-9)
The flux computed by SOILINER was then used to compiute Gardner's solution
(4-8) using the computer program listed in Appendix B.
Figure 4-1 shows the agreement between SOILINER simulations with a
regular grid (Test 2A) and a graded grid (Test 2C) and Gardner's analytic
solution.
1-35
-------
09
I
CD
(VI.
z: •
N
CD
II).
TEST Efl
TEST EC
flNflLYTIC
-B-
.0
.5
1.0
1.5
PF
2.0
2.5
3.0
Figure 4-1. Comparisons of solutions for steady state vertical flow
upward from a water table: solid curve is analytical
solution of Gardner (1958); symbols are SOILINER numerical
solutions at every other node. Test 2A (D) is a regular
spaced grid with 50 nodes, Test 2C (+) is a variable spaced
grid with 46 nodes. (pF = log (-^), ^ in cm.)
1-36
-------
4.3 Transient Infiltration
Philip [1958, 1969] developed a quasi-analytical solution to transient
infiltration into unsaturated soil and applied this technique to simulate
infiltration into the Yolo light clay with a porosity of 0.495 and a saturated
hydraulic conductivity of 1.23 x 10"-* cm/s. The functional relationships
for the unsaturated hydraulic conductivity and moisture characteristic were
taken from Haverkamp et al. [1977] and are shown in Milly [1982]. The initial
boundary conditions for this problem are
ijj = -600 cm all z, t <_ 0
4> = 25 cm z = 0, t > 0
* = -600 cm z = - °°, t > 0
The soil is initially at a constant moisture content and capillary pressure
over its entire thickness. For the SOILINER simulations, the initial time
step was 0.1 sec and subsequent time steps were automatically calculated to
keep the maximum pressure change between time steps at any node to about 10 cm.
Figure 4-2 shows a comparison of SOILINER using a regular grid with
Philip's quasi-analytic solution. The agreement for this highly nonlinear
problem is very reasonable. Figure 4-3 shows the similar accuracy of SOILINER
with a graded grid. The graded grid has fewer nodes, but because of smaller
node spacings at the top of the column, the initial moisture profile is closer
to Philip's results, demonstrating the value of variable grid spacing.
SOILINER is shown to accurately simulate the flow of moisture in a
vertical unsaturated column using both regular and graded grids. The
infiltration test is very similar to the application of SOILINER in liner
design. Successful completion of these verification runs using other computer
models would similarily indicate their utility for soil liner design.
4.4 Particle Tracking
To verify the particle tracking algorithm, a fully saturated clay liner
was simulated using SOILINER for comparison with a calculated breakthrough
time based on Darcy's Law (equation 2-2). The following boundary and initial
conditions were imposed on a 90 cm clay liner having a conductivity of
1.0x10"^ cm/sec and porosity of 0.495: (1) impoundment depth of 100.0 cm,
(2) pressure at the liner base of 0.0 cm, and (3) a steady state pressure
distribution within the liner at time t = 0. As required for steady state
conditions, the change in flux over space is zero (i.e., flux everywhere in
the liner is constant - see equation 2-lb). Since the entire clay liner is
saturated, moisture content is also constant. Thus, regardless of the
particle location in the liner, its velocity will remain unchanged.
Figure 4-4 reveals this constant velocity up to the point of breakthrough at
year 6.69.
1-37
-------
MOISTURE COMTEMT
N
25 .38 .35 .48 .45 .50
Figure 4-2.
Comparison of Philip's quasi-analytic solution (a)
and SOILINER regular grid solution (solid line) for
infiltration into Yolo light clay under ponding.
Curves are at 103, 104, 105, and 2xl05 sec.
1-38
-------
MOISTURE CONTEMT
.28 .25 .30 .35 .48 .45
I i i | I
.58
CD
6)
I
CD
l
z:
u
N
CD
tn_
CD
*-.
I
CD
U)
Figure 4-3. Comparison of Philip's quasi-analytic solution (o) and
SOILINER graded grid (46 nodes) solution (solid line) for
infiltration into Yolo light clay under ponding. Curves are
at 103, 104, 105, and 2xl05
sec.
1-39
-------
PARTICLE TRACKING VERIFICATION
.p-
o
E
o
Time (yrs)
Figure 4-4. Change in particle depth with time across a fully saturated
clay liner under steady state conditions.
-------
For comparison, an estimated breakthrough time under the same initial and
boundary conditions (across the entire liner thickness) can be obtained from
equations 3-23 and 3-24:
q = 1.00xlO~7 cm/sec *
U-10Q.O
-90.0
+ 1
= 2.11xlO~7 cm/sec
v = 2.11xlO~7 cm/sec * 0.495
= 4.26xlO~7 cm/sec
t = 90 cm * 4.26xlO~7 cm/sec
= 2.11xl08 sec 3.1536xl07 sec/yr
= 6.69 years
This calculation compares accurately with the breakthrough time of 6.69 years
obtained from S01LINER.
1-41
-------
PART II
SOILINER USER'S MANUAL
II-l
-------
1. INTRODUCTION
1.1 History of the SOILINER Computer Code
The SOILINER computer model was originally developed at GCA under
contract to EPA's Office of Solid Waste, Land Disposal Branch. The first
version was completed in September, 1983. In April, 1984 the SOILINER report
was released by EPA as a Technical Resource Document (EPA/530-SW-84-001; Goode
and Smith, 1984) for public comment. This document emphasized the theoretical
development of the finite-difference technique utilized in SOILINER and the
ability to accurately simulate an infiltration event.
In November 1984 GCA investigated the times to breakthrough and steady
state, with associated fluxes, for a variety of clay liner scenarios using
SOILINER (Johnson and Wood, 1984). At that time, SOILINER was not
specifically coded for the determinations of breakthrough and steady state.
However, data obtained from the calculated moisture contents provided by
SOILINER served as the basis for determining these times. The results of that
study indicated that a moisture-based criterion for establishing breakthrough
was not appropriate (see Section 3.4 of Part I) and a less subjective method
should be developed.
The fundamental basis for simulating an infiltration event remains
unchanged in the newest version of SOILINER. However, the original model has
beevi modified to determine the times at which steady state and breakthrough
occur. Breakthrough is determined from a particle tracking algorithm which is
discussed further in Section 2.
1.2 Scope of the Manual
The remainder of this document is devoted to the development of a
step-by-step procedure for the application of SOILINER. Section 2 discusses
in detail the features and limitations of SOILINER such that it becomes
apparent which flow scenarios can be accurately simulated. Section 3 provides
a general description of the data necessary as input to the model. This
section also discusses the development of characteristic curves for moisture
content and unsaturated hydraulic conductivity, both a function of matric
potential. It is important to note that these two relationships form the
basis of an accurate simulation. Although a number of characteristic curves,
covering a wide range of soil types, have already been coded into the model,
site specific data can be easily incorporated into SOILINER.
Section 4 presents the findings from a thorough testing of the SOILINER
model for a variety of solution strategies. Testing includes the
investigation of modifications to the original code and new programming
features. A number of simulations have been made, from which specific
guidelines were established describing appropriate values for the model
parameters for a variety of clay-liner scenarios.
II-2
-------
1.3 Computer Requirements
SOILINER was originally written in FORTRAN IV for an IBM 3033 mainframe.
The current single-precision version is dimensioned for a maximum of 200 node
points. The average run-time for a 100_+ node flow domain over a considerable
time period (3 to 5 years) takes between 30 and 60 seconds of CPU-time.
A modified version of SOIUNER has been run successfully on an IBM-PC
(256 K), and requires approximately 128 K-bytes of RAM. Even with an
Intel 8087 math co-processor, average run-times are considerably longer.
Both versions require 4 input files and 5 output files. Input files
include the following:
(^ Grid File - configures node coordinates.
(2) Properties File - configures element soil properties by layer,
including the saturated properties of the soil from which the
characteristic curves were developed.
(3) Initial Conditions file - configures initial fy distribution.
(4) Control File - contains parameters which define the chosen solution
strategy including the number of nodes, time stepping parameters,
and boundary conditions.
Output files include:
(1) General Output - this file includes an echo of input data and
pertinent output by time step including moisture, pressure, and
conductivity results.
(2) Matric Potential, and (3) Moisture Content - these two files repeat
data from the General Output File but are reproduced separately for
graphical analysis.
(4) Flux Output - this file contains element flux data by time step, and
can be used to demonstrate flux divergence (note that negative flux
values indicate downward movement, out of the liner whereas positive
values indicate movement up, into the liner).
(5) Velocity Output - this file contains data relevant to particle
tracking and can be used to represent the dynamics of infiltration
and the time to breakthrough.
The use of separate files facilitates input data manipulation and the
graphical analysis of output. Individual parameters and their format are
discussed further in Appendices C and D.
II-3
-------
2. MODELING CAPABILITIES
2.1 Model Assumptions
SOILINER was specifically developed to simulate vertical, unsaturated
flow beneath ponded liquid. The assumptions implicit in the governing
equation (2-1) are:
(1) the air phase is always at a spatially constant atmospheric pressure;
(2) the fluid has a constant density and does not freeze;
(3) flow is not affected by temperature gradients or solute
concentration gradients; and
(4) the medium does not deform.
The last assumption must be carefully considered, especially for liners having
a high percentage of clay which may be subject to leachate interaction. Under
extreme conditions, dissolution and piping (the actual loss of clay particles
to underlying strata) may occur due to the chemical nature of an impounded
liquid. Other constituents have been shown to alter the structure of clay
barriers. In any case, liner deformation invalidates the characteristic curve
for the soil under consideration. However, most hazardous waste leachates are
water-based. While strongly concentrated leachates may increase permeability,
no examples were found in the literature where very dilute aqueous solutions
caused substantial permeability increases (Anderson and Jones, 1983).
One final consideration is specific to the numerical representation of
the governing equation (2-1). Most of the characteristic curves available in
the SOIUNER model were developed utilizing data obtained from the desorption
cycle. Consequently, SOILINER does not account for hysterisis.
2.2 Solution Strategies
SOILINER is capable of running transient and/or steady state solutions
depending on the combination of parameters established in the control file.
If a steady state solution only is desired, the ISTEDY integer flag must be
set equal to 1. For a transient simulation ISTEDY is set equal to 0, and a
MAXNT is established such that all desired output times are reached for the
given combination of time stepping parameters (see Section 4). If the
transient option is chosen, the steady state algorithm is first solved as a
means of determining the time at which steady state is achieved during
transient simulations. Additional parameters for both the steady state and
transient solution strategies are covered in Section 4 and Appendices C and D.
II-4
-------
2.3 Definition of Clay Liner Domain
The size and content of the clay liner domain depends on the limits
imposed on the system such that boundary conditions can be specified. The
upper boundary should always be specified at the clay-liner upper surface and
will generally have \\> set equal to the impoundment depth. Ideally the lower
boundary would be established at the water table where i|> = 0. This domain may
include a large portion of the underlying site soil, and would require
discretization of both the liner and native soil. However, it is possible to
model only the clay liner if the lower boundary condition can be reasonably
estimated. When a sharp discontinuity in hydraulic conductivity exists
between the liner and native soil (i.e. 1.0 x 10~' to 1.0 x 10 cm/sec),
the lower boundary can be set to reflect the height above a water table. For
comparison, Figure 2-1 depicts the changes in if over the period of
infiltration for a 90 cm clay liner (K - 1.0 x 10~7 cm/sec) underlain by
300 cm of an unsaturated, sandy soil (K = 9.4 x 10"-* cm/sec). Figure 2-2
depicts the ^ distributions for a similar domain where the only difference is
a decreased conductivity for the underlying, native soil
(K = 1.0 x 10~5 cm/sec). As shown in Figure 2-1, the initial ifi of -300 cm
at the clay/sand interface increased only 9 cm to -291 cm, whereas the ^
increased from -300 cm to -137 cm if the sand were replaced with a silt (see
Figure 2-2).
SOILINER is also capable of modeling layered flow domains. To
demonstrate this ability, a two liner system was simulated using SOILINER.
The hypothetical liner was constructed as follows:
Layer 1: 61 cm (21) clay; K = 1.0xlO~7 cm/sec, n = 0.495
Layer 2: 30 cm (I1) sand; K = 9.44xlO~3 cm/sec, n = 0.287
Layer 3: 91 cm (3' ) clay; K = l.OxlO"7 cm/sec, n = 0.495
Native Soil: 300 cm to water table (PSI = 0.0); K = 9.44xlO~7 cm/sec,
n = 0.287
where the impoundment depth was set equal to 100 cm. To simulate worst case
conditions, a leachate collection system was not included in the upper sand
layer.
The change in pressure with time shown in Figure 2-3 reveals the dynamic
interaction between layers in this liner configuration. Figure 2-4 depicts
the corresponding changes in moisture content. Note that a drop in if*
II-5
-------
E
o
c
0)
~o
CL
O
100
50 -
0
-50 -
-100 -
c -150
-200 -
-250 -
-300
100
200
Depth (cm)
300
400
Figure 2-1. Changes in pressure distribution for a clay liner underlain
by a high conductivity native soil.
-------
E
o
c
0)
4->
o
Q.
O
100
50 -
0
-50 -
-10O
r -150
-200 -i
-250 -
-30O
Q INITIAL CONDITIONS
1OO
200
Depth (cm)
300
400
Figure 2-2. Changes in pressure distribution for a clay liner underlain
by a low conductivity native soil.
-------
OC
E
o
c
-------
C
o
o
0)
D
-*->
(0
'o
POINT OF CLAY
SATURATION
(0 = n =0.495
POINT OF SAND
SATURATION
(0= n =0.287)
DINITIAL CONDITIONS
4- YEAR
OYEAR
A YEAR
X YEAR
V YEAR
0.5
.0
.0
3.0
3.8
(STEADY STATE)
200
Depth (cm)
Changes in moisture content with time corresponding
to pressure curves shown in Figure 2-3.
400
-------
indicates moisture loss, whereas an increase in \\> indicates moisture gain.
Many salient features of the infiltration event can be drawn from these
graphics:
(1) The upper portion of the first sand layer loses moisture at a rate
greater than replenishment from the overlying clay (see year 0.5
where ij; in the sand drops below initial conditions); this net loss
is offset by moisture gain in the clay immediately below the sand.
(2) For the first three time steps moisture is lost from the bottom clay
liner base, as indicated by the ^ drop which propagates into the
liner.
(3) Overall, pressure becomes less negative as pores reach their steady
state moisture content; note that pores initially drained at the
liner base are rewetted from infiltrating liquid sometime after
year 2.
(4) Although ip increases with depth in both sand layers, the total head
OJ>+ z) always decreases with depth resulting in downward flux.
(5) At steady state, the entire two-liner system (181 cm) is saturated
to a depth of 164 cm, at which point ty values drop below 0.0 with
increasing depth.
For this liner configuration with an initial pressure of -150 cm, steady
state was achieved in 3.8 years and particle breakthrough occurred at
year 11.5.
2.4 Determining Breakthrough
SOILINER utilizes a particle tracking algorithm to determine the time
required for a particle to be advected across a specified liner thickness.
The effects of dispersion are not accounted for, and thus, breakthrough
represents the average linear velocity of migrating contaminants. The liner
thickness of concern is established by setting NCLAYN equal to the bottom node
of the liner. The particle is initially positioned at the liner surface.
Subsequent movement is then determined from velocity data generated by
SOILINER. Figure 2-5 depicts the dynamics of particle movement across a 60 cm
clay liner. The initial velocities are high as shown by the steeply sloping
curve near the liner surface (i.e., Az great with respect to At). With time,
velocities decrease as the wetting front penetrates the liner.
2.5 Variable Initial and Boundary Conditions
Use of the numerical technique not only enables one to vary the soil
properties, but also allows variable initial conditions to be specified by the
11-10
-------
E
u
Q.
- c
0.0
-10.0 -
-20.0 -H
-30.0 -
-40.0 -
-50.0 -
60.0 -
-70.0
Particle Location (61cm Clay Liner)
Time (yrs.)
Figure 2-5. Changes in particle depth depicting the dynamic characteristics of the infiltration
event, specifically decreasing wetting front velocity with respect to time.
-------
user. It may be reasonable to assume a constant pressure for a carefully
constructed, homogeneous clay liner in which the moisture content for each
lift is monitored and adjusted. However, multilayered liners are likely to
exhibit widely varying initial pressures and associated moisture contents.
Also, the pressure distribution in native soils prior to liner construction
will vary, generally becoming more negative with increasing height above the
water table (see Section 2, Part I). By breaking the continuous flow domain
into a limited number of discrete node points for the finite-difference
technique, it is possible to specify different initial pressures at various
depths in the liner domain. This inherent ability also facilitates
investigation of the effects on infiltration due to time-varying boundary
conditions.
At any specified point in time (t^) it is possible to predict the i|>
pressure distribution for a given set of initial and boundary conditions. At
that time it may also be desirable to determine the effect of an increase or
decrease in impoundment depth (i.e. upper boundary condition). To achieve
this objective with SOILINER, the user would simply employ the existing,
calculated pressure distribution at time t^ as initial conditions for
subsequent simulation under the new impoundment depth. Once the new boundary
conditions are specified, a desired set of special output times can be
obtained which will reflect the changes imposed on the system.
Figure 2-6 depicts the initial fy distribution and four special output
times for a clay liner simulation where the impoundment depth was changed from
100 cm to 200 cm prior to particle breakthrough. For the first three output
times (up to year 1.5) the impoundment depth was constant at 100 cm. By
year 1.5 steady state was achieved and the particle had penetrated the liner
to a depth of 36 cm. At this time the impoundment depth was increased to
200 cm. Within 2 months, a new steady state pressure distribution was reached
(year 1.7).
The impact of increasing impoundment depth on particle movement is
depicted in Figure 2-7. If the head (H) had remained constant at 100 cm, as
shown by the dotted line, breakthrough would have occurred at approximately
year 5. However, the change in head (AH) at year 1.5, increased the particle
velocity due to an increased gradient, with particle breakthrough occurring
one year earlier.
II-12
-------
20O
VARIABLE HEAD IN IMPOUNDMENT
1OO
QINITIAL CONDITIONS
+ MONTH 3.0
O MONTH 9.0
A YEAR 1.5
V YEAR 1.7
E
o
c
0>
-f-'
o
Q.
O
O
-10O -
-2OO -
-300
n i 11 u i n n n n n
0
20
Figure 2-6.
40
60
80
100
120
Liner Depth (cm)
Change in steady state pressure distribution due
to increased impoundment depth at year 1.5.
-------
VARIABLE HEAD IN IMPOUNDMENT
E
o
JC
-**>
CL
Q
Time (yrs)
Figure 2-7. Effect of increased impoundment depth on
the time to breakthrough.
-------
3. DATA REQUIREMENTS
Data necessary to run SOIUNER are related to either (1) the solution
strategy, or (2) liner properties. Most of the parameters found in SOILINER
are required solely for the development of a solution strategy (i.e. time
stepping, grid discretization, etc.) to simulate a given set of conditions.
Although very little field data is required, the method in which soil
properties are obtained (and their accuracy) is critical and serves as the
fundamental basis for the SOILINER model. Methods of determining soil
moisture and conductivity as a function of matric potential (i.e. negative
pressure) will be described in detail after a general discussion of the
required input.
3.1 Description of Input Data
As mentioned previously, there are four input data files. Following is a
listing, by file, of the input data.
Grid File:
NUM - number of consecutive nodes with separation DELTA.
DELTA - separation between adjacent nodes (cm).
Properties File:
NUM - number of elements in soil layer.
ISOIL - one-dimensional array holding soil-type code numbers
ranging from 0 to the number of soils for which
characteristic curves have been coded into the model (see
Table 3-1, p. 11-25).
SATK - saturated conductivity for the associated soil type.
FOR - porosity of the associated soil type.
PSICRT - critical pressure (see Section 3.3, Part II).
Initial Conditions File:
NUM - number of consecutive nodes with initial pressure, PSI.
PSI - initial pressure, PSI (cm).
Control File:
TITLE - character variable used to describe each particular run.
11-15
-------
LPRINT - Boolean flag for output format; X = long, F = short.
LSTEDY - Boolean flag; T = steady state solution, F = transient.
ENDTIM - desired total period of simulation (yrs).
DT - initial time step size (sees).
CHPARM - change parameter utilized to increment DT between successive
time steps.
DTMAX - maximum allowable time step (sees).
MAXNT - maximum number of time steps.
ALPHA - time weighting parameter set between 0.0 and 1.0 for
transient solutions.
ERRMAX - criteria established for the maximum change in ^ between
successive iterations at any given time step (cm).
MAXIT - maximum allowable number of iterations at any given time
step.
NOUT - number of special output times.
TOUT - one-dimensional array holding NOUT number of special output
times (NOTE: if NOUT <0, TOUT must not be present; if
NOUT = -1, all time steps are printed).
NUMNP - total number of nodes.
NCLAYN - number of liner nodes (utilized in determining breakthrough),
SRPARM - successive relaxation parameter ranging from 0.0 to 1.0.
H - depth of liquid impoundment (upper ^ boundary condition).
PSIBOT - pressure at the base of the flow domain (lower ^ boundary
condition).
Review of the necessary input data reveals that the only field- or
lab-oriented measurements are the initial pressures and the soil properties
SATK and POR. What is not readily apparent, nor required as direct input from
11-16
-------
the user, are the characteristic curves for the specified soil types. The
parameters H, PSIBOT, initial PSI, SATK, and POR are readily available (either
measured or estimated) in comparison to the effort required to develop a
characteristic curve. Subsequent reference to parameters listed in the above
files will utilize their FORTRAN variable names when appropriate. Appendix C
provides a more complete list of variables found in S01LINER.
3.2 Measurement of the Soil Moisture Characteristic Curves
In measuring the fundamental relationship between soil moisture and
suction (raatric potential x specific weight of water) it is better, in
principle, to take in-situ measurements rather than measurements of disturbed
samples (e.g., dried, screened, and artificially packed sample) in the
laboratory. At low suction (0-1 bars) the soil moisture retention is strongly
influenced by soil structure and pore size distribution (Hillel, 1980). The
most preferable method in measuring the soil moisture characteristic is by
making simultaneous measurements of wetness, with a neutron moisture meter,
and suction, with a tensioraeter. Unfortunately, this approach has often been
frustrated by soil heterogeneity and by uncertainties over hysteretic
phenomena as they occur in the field. The following paragraphs briefly
describe these instruments for use in measuring soil moisture and suction in
the field. The reader is referred to Hillel (1980) for more detail.
The tensiometer is the most widely accepted practical device for
measuring matric potential in the field. As shown in Figure 3-1, the
tensiometer consists of a porous cup connected through a tube to a manometer,
where all parts are filled with water. The porous cup, which is generally
made of ceramic material, is placed in the soil where the suction measurement
is to be made. When initially placed in the soil, the water contained inside
the tensiometer is generally at atmospheric pressure. The soil water
surrounding the cup is most likely below atmospheric pressure. This causes a
suction which draws a certain amount of water from the tensiometer, thus
causing a drop in its hydrostatic pressure. This pressure is registered on
the manometer, which can be a simple water or mercury filled U tube, a vacuum
gauge, or an electrical transducer.
A tensioraeter is placed in the soil for a long period of time to follow
changes in the pressure of soil water. As moisture is depleted from the soil
by evapotranspiration or drainage, or replenished by percolation or
irrigation, corresponding readings on the tensiometer can be obtained.
During measurements of suction, simultaneous readings of soil moisture
content must be obtained at the same locations in order to relate soil
moisture with suction. An efficient and reliable instrument for measuring
soil moisture is the neutron moisture meter. As shown in Figure 3-2, this
instrument consists of a probe which is lowered into the soil, attached to a
scalar or ratemeter. The probe contains a source of fast neutrons (obtain an
average speed of 1,600 km/sec) and a detector of slow neutrons. The scalar
monitors the flux of slow neutrons scattered by the soil.
11-17
-------
In general, the probe emits fast neutrons into the surrounding soil where
they encounter and collide with the atomic nuclei of hydrogen. When they
collide, the fast neutrons lose some of their kenetic energy, thus becoming
slow neutrons. Neutrons that have been slowed to a critical speed are said to
be thennalized. As the thermalized neutrons repeatedly collide, and move
about randomly, a number of them return to the probe where they are counted by
the slow neutron detector. The amount of thermalized neutrons is proportional
to the concentration of soil moisture.
3.3 Functional Relationship Formula
To formulate the two part function
4> = *SW"4 for 0 < W < W£ (3-1)
and
ty = -m(W-n)(W-l) for Wi < W < 1 (3-2)
where m and n are:
m
(i-w.)2
(3-3)
n - 2W. - —- - 1
x mW.
i
in which the parameters ijj, i|>8, W, b, V^, and 4^ are shown in Figure 3-3,
soil moisture content, 8, and matric potential, 4>, must first be obtained
using field techniques described in the previous section or from the
literature (note that ^i and Qs are PSICRT and FOR, respectively).
To calculate soil wetness, W, the unsaturated moisture content, 9, must
be divided by the saturated moisture content, 9S or in this study, the total
porosity. The saturated moisture content must be determined by first oven
drying an undisturbed sample from the field, and then weighing it. It is then
saturated with some liquid and weighed again. Finally, the saturated sample
is immersed in the same liquid, and the weight of the displaced liquid is
noted. The weight of the liquid required to saturate the sample divided by
the weight of the liquid displaced is the porosity as a decimal.
Both the matric potential at saturation, i^s, and the exponent b are
empirical and are estimated using a power curve regression analysis performed
on the log transforms of the obtained data (4 and W, where W = 8/6g). The
value 4^3 is the log of the y-intercept, and b is the slope of the line.
To determine the values of the parameters m and n, the location of the
inflection point, (W^, 4>j_) in Figure 3-3 must be determined. This is
conducted by first plotting the soil moisture matric potential data and
11-18
-------
Vacuum
gaui
Soil surface
Openi ng
tC f : !' *.1
water
Air
trap
Connect -
ing tube
Porous
ceramic
cup
U
a>
Depth - d
Figure 3-1. Illustration of the essential parts of a tensiometer
(from Hillel, 1980) .
:AELE
FACE
OAFJlOACTUE SOURCL
•-ACTIVE ZONE
ACCESS TUBE
Figure 3-2. Illustration of the essential components of a portable
neutron soil-moisture meter (after Hillel, 1980).
11-19
-------
b = slope of log —log plot of
versus W for 0
-------
observing the point where d^/dW changes from an increasing to a decreasing
function as W decreases. It should be noted that equation 3-2 represents the
W versus ty curve where W^< W <1 only if m > 0, which requires that W^>
b/b(b+l). Once the inflection point has been located, the values i]^, W^,
and b are used to calculate m and n using equations 3-3 and 3-4. In turn, tyai
b, m, and n are inserted into equations 3-1 and 3-2 to obtain the two
part function described by Clapp and Hornberger (1978) relating soil moisture
and matric potential.
Rearranging equation 3-2 using the quadradic formula gives
J~
W =
which is a convenient arrangement for solving W in the FORTRAN computer code,
given m, n, and fy are known.
The power curve function
k = W2b+3 (3-6)
is used in the SOILINER model to relate soil moisture, 6 (W = 6/6s), and the
unsaturated hydraulic conductivity, K, where
k = K/KS
in which k is the relative hydraulic conductivity, Ks is the saturated
hydraulic conductivity (SATK), and b is the same value as the parameter b in
equation 3-1. To formulate this function for a specific soil, Ks must be
determined.
Measurements of a saturated hydraulic conductivity can be made in the
laboratory using falling or constant head perraeameters, or in the field by
means of a bail or slug test. In using permeameters, a sample of the soil is
subjected to water under a known or falling head, and the flow through the
sample in a known time is measured. Such tests have limited value due to the
difficulty of placing samples representative of their natural state in the
permeameter. Also, flow in solution cavities cannot be duplicated in a
permeameter.
Bail or slug tests are carried out in the field using a single
piezometer. Both tests are conducted by causing an instantaneous change in
the water level within the piezometer. This change is achieved by either
adding a known volume of water into the piezometer (slug test) or, removal of
a known volume of water from the piezometer (bail test). The recovery of the
water level is then observed.
A number of methods have been derived for determining the saturated
hydraulic conductivity from the recovery data. The simplest interpretation of
the recovery data is that of Hvorslev (1951). This method is for a point
11-21
-------
piezometer. For bail or slug tests run in piezometers that are open over the
entire thickness of an aquifer, the methods of Cooper et al. (1967) and
Papadopoulos et al. (1973) should be used. The main limitation of slug or
bail tests is that the well point or screen should not be clogged or
corroded. In addition, it may be difficult finding an in-situ, saturated soil
that has the same characteristics as the liner soil being modeled.
Once the values of ^s, Ks, 0S, b, W^, and 4^ have been
determined, equations 3-1, 3-5, and 3-6 are inserted into subroutine SPROP in
the SOILINER model (Appendix G). A soil type number is assigned to the
particular soil so that it can be distinguished from other soil types which
currently exist in the model.
3.4 Characteristic Curves Available in SOILINER
The SOILINER model uses functional expressions for determining
unsaturated hydraulic conductivity, moisture content, and specific moisture
capacity (which is the derivative of moisture content with respect to matric
potential). For each element within the soil column, the soil type code (1,
2, 3, etc.) must be selected and used as input in the Soil Properties File.
This allows stratification of different soil types within a column (i.e.,
clay, silt, sand, clay). The FORTRAN listing shows that 13 soil types are
used, but the number of soil types is limited only by the number of functional
relationships inserted into the subroutine SPROP. Table 3-1 shows the soil
types used in the model and their associated values for saturated hydraulic
conductivity, saturated moisture content (equal to the porosity in this
study), and saturated matric potential. Figures 3-4a and 3-4b show the soil
moisture characteristic curves for soil types 1 through 11 listed in Table 3-1,
The soil types numbered 1 through 11, and their associated empirical
equations in Appendix G, are from Clapp and Hornberger (1978). Of these
soils, the silty clay has the lowest saturated conductivity value of
1 x 10"^ cm/8. Typically, a clay liner consists of clay with a saturated
hydraulic conductivity of lxlO~7 cm/a (40 CFR Part 264). Based on an
extensive literature search, GCA was unable to locate any characteristic data
representative of a clay with such a low saturated conductivity. Therefore, a
hypothetical clay soil was derived by using the Yolo light clay utilized by
Haverkamp et al. (1977) with a modified saturated conductivity of
1 x 10"^ cm/s (soil type 12). The empirical equations derived by Haverkamp
et al. (1977) to relate matric potential with soil moisture and unsaturated
hydraulic conductivity are also used in the SOILINEK model (Appendix G) to
represent this hypothetical clay. Finally, soil type 0 represents a high
conductivity sand also developed by Havercamp et. al. (1977).
11-22
-------
C
0>
-»-•
C
O
u> • —
^
0)
o:
0.9 -
PF=(LOG(-PSI))
Figure 3-4a.
Characteristic moisture curves for soil types ranging from
sand to sandy clay loam.
-------
I
fo
c
(0
-t->
c
0
(J
Q)
0)
>
0
PF=(LOG(-PSI))
Figure 3-4b.
Characteristic moisture curves for soil types ranging
from silty clay loam to a compacted clay.
-------
TABLE 3-1. REPRESENTATIVE VALUES OF HYDRAULIC PARAMETERS FOR
THE SOIL TEXTURE CLASSES CODED INTO SOILINER
Soil code
number
0
1
2
3
4
5
6
7
8
9
10
11
12
Soil
texture
Sand
Sand
Loamy sand
Sandy loam
Silt loam
Loam
Sandy clay
Silty clay
Clay loam
Sandy clay
Silty clay
Clay
Heavy clay
Mean clay
fraction
—
0.03
0.06
0.09
0.14
0.19
loam 0.28
loam 0.34
0.34
0.43
0.49
0.63
**
SATK (K )
(cm/sec)
9.44xlO~3
1.760xlO-2
1.563xlO~2
3.466xlO~3
7.200xlO~4
6.950xlO~4
6.300xlO~4
1.700xlO-4
2.450xlO-4
2.166xlO"4
1.033xlO~4
1.283xlO~4
l.OOOxlO"7
POR (9 )
s
0.287
0.395
0.410
0.435
0.485
0.451
0.420
0.477
0.476
0.426
0.492
0.482
0.495
PSICRT (4*. )
(cm/sec;
-1.00
-16.96
-12.97
-32.08
-122.28
-74.92
-54.12
-67.94
-128.19
-36.42
-116.63
-104.78
-1.00
* From Li et al. (1976)
** Hypothetical clay liner soil from Havercamp et. al. (1977)
11-25
-------
4. PROCEDURE FOR APPLYING SOILINER
Upon completion of the initial data gathering phase, it should be
possible to synthesize a conceptual model of the liner. A conceptual model is
an abstract description of the flow domain including the: (1) liner
configuration, (2) nature and position of the boundary conditions, and
(3) initial pressure distribution. Once the conceptual model has been
established for a desired simulation, SOILINER can be applied.
4.1 Model Set-up
Following is a brief, step-by-step procedure in logical sequence to
establish a SOILINER simulation:
1. Locate boundaries of the flow domain based on the liner thickness
and depth to the water table. Boundaries must be positioned such
that values of ty (H and PSIBOT) can be specified.
2. Design a finite difference grid (discretization). Appendix D - Grid
File provides an example discretization scheme.
3. Choose an appropriate characteristic curve for each soil type and
assign the corresponding saturated properties (see Section 3.3) to
each element of the flow domain as shown in Appendix D - Properties
File.
4. Construct the initial pressure distribution based on the design
specifications of the liner material during installation (i.e.,
optimum moisture content and the corresponding matric potential) and
the existing field conditions for the native soil as described in
Appendix D - Initial Conditions File. Typically, the moisture
profile in the native soil is assumed to be in static equilibrium
above a fixed water table which is deemed to represent an average
depth below the liner base (see Section 2, Part I).
5. Develop a solution strategy as shown in Appendix D - Control File.
Steps 1, 3, and 4 are fairly straightforward and have been discussed in
detail previously. More attention will be given to methods of discretization
(Section 4.2) and the development of solution strategies (Sections 4.3 and
4.4). The discussions on steady state and transient solution strategies are
based on an in-depth sensitivity analysis of parameters within SOILINER.
4.2 Designing a Finite Difference Grid
Once the dimensions of a conceptual model have been established, it is
necessary to discretize the flow domain. Discretization results in a finite
number of node points which constitute a one-dimensional grid. The governing
equation employed by SOILINER is solved in terms of these node points (see
Section 3, Part I). Thus model accuracy is influenced by the discretization
11-26
-------
scheme, and the following must be considered during grid development:
(1) locations of abrupt changes in the initial PSI distribution, usually
occurring at the liner surface, (2) sharp discontinuities in soil properties,
and (3) model efficiency - although a large number of nodes will generally
increase simulation accuracy, computer storage requirements and associated run
times will increase.
The simplest grid divides the liner and site soil into equally sized
subdoraains. These elements may be fairly large if there are only small
variations of moisture content in the profile, such as occurs under conditions
of equilibrium. However, if large gradients of moisture (e.g., a wetting
front) occur, the size of the elements should be smaller. In order to
economize on the total number of elements, it is important in such cases to
use a variably-sized set of elements in order to concentrate nodes in the
region of greatest variations, as discussed below.
For the transient infiltration problem, the initial distribution of
matric potential varies gradually everywhere, except at the top of the liner.
As infiltration proceeds, this sharp front moves into the soil and is
gradually smoothed as it moves down the soil column. By the time the moisture
front reaches the bottom of the liner, it is sufficiently dispersed to relax
the subdomain size requirement. This variation suggests a grid with small
subdomains at the top of the liner and a gradation of subdomain size moving
down the column. Figure 4-1 shows two graded grids of this general design.
In the first, the grid consists of blocks of segments which become larger in
steps. The second sample has subdomains which vary in size continuously from
the liner top to the water table. In this second grid, the bottom subdomain
is about 10 times as large as the top one.
It is difficult to specify a general criterion for the size of grid
spacing. The subdomain size and gradation are usually determined during the
initial simulations. Many models are very sensitive to this discretization
and care should be taken. A good test of the effect of subdomain size on the
solution is to compare the change in simulated matric potential obtained using
smaller and smaller grid spacings. When the change between two simulations is
unaffected by grid spacing, the larger size can probably be used.
4-3 Steady State Solution Strategy
Regardless of the solution strategy chosen, SOILINER will always
calculate the steady state PSI distribution. For transient simulations, the
steady state solution serves as a means of determining the time to steady
state during an infiltration event. However, it is possible to request only
the steady state solution (by setting LSTEDY = T) for the purpose of first
developing this strategy. This feature avoids lengthy run times associated
with transient simulations before properly developing a steady state solution.
In comparison with the transient formulation, there are relatively few
parameters in the Control File which are required for a steady state
simulation - SRPARM, MAXIT, and ERRMAX. Of these three, only the successive
relaxation parameter (SRPARM) is specific to the steady state algorithm.
11-27
-------
BLOCK GRADED
CONTINUOUS GRADATION
T)
BLOCK 1
AZ =
BLOCK 2
AZ=2AZ1
BLOCK 3
AZ=5AZ1
SUBDOMAIN 1
SUBDOMAIN NUMEL
AZ;
AZ
AZNUMEL
NUMEL = number of subdomoins
or elements
Figure 4-1. Two vertical grids with variable subdomain sizes.
11-28
-------
As discussed in Part I, the steady state solution is achieved
iteratively. For each iteration the calculated DSPI at each node is forced to
satisfy the specified boundary conditions based on the current soil
properties. The pressure distribution is then updated as follows:
PSI = PSI + (SRPARM * DPSI) (4-1)
Soil properties are adjusted to reflect the new PSI distribution and DPSI is
again calculated. The iterative procedure continues until the change in soil
properties no longer significantly affects the calculated PSI distribution
(i.e., DPSI < ERRMAX).
For layered systems, values of SRPARM > 1.0 lead to model instability.
Figure 4-2 depicts the affect of SRPARM = 1.0 on the change in PSI at nodes
25, 55, and 85 from a 500 era liner configuration (180 cm clay, 320 cm sand).
Without relaxation (SRPARM = 1.0), the magnitude of the calculated DPSI
between iterations is far greater than a reasonable ERRMAX at each node.
Similar oscillations occur at all other nodes. To reduce numerical overshoot
and the propagation of subsequent oscillation, it is necessary to specify a
value of SRPARM that is less than 1.0.
When SRPARM < 1.0, only a percentage of the calculated DPSI is used to
resolve PSI during the iterative procedure. Thus, the magnitude of initial
overshoot is reduced and subsequent oscillations may be sufficiently damped to
allow convergence. Figure 4-3 reveals the affect of SRPARM< 1.0 on
convergence at node 85. As shown, SRPARM = 0.80 was not sufficient to achieve
convergence. However, values of 0.6, 0.4, and 0.2, all converged on the same
steady state PSI distribution. Termination of the iterative procedure at
node 85 for these three values indicates that all nodes in the flow domain
also met the error criterion.
Although values for SRPARM of 0.2, 0.4, and 0.6, all lead to convergence,
the number of iterations did vary (27, 16, and 24, respectively). Tables 4-1
and 4-2 provide a comparison on the efficiency of various SRPARM scenarios for
a 90 cm and 180 cm clay liner configuration, respectively. For both
scenarios, MAXIT was set at 100 and ERRMAX was established at 0.01 cm.
Results indicate that: (1) SRPARM generally works best at values close to
0.5, (2) decreasing head values require more iterations, (3) increasing MAXIT
for values of SRPARM > 0.5 does not guarantee convergence (these scenarios
frequently diverge if SRPARM > 0.5), and (4) as the number of nodes increases,
more iterations are required before ERRMAX is satisfied simultaneously at all
nodes. In general, for all cases where SRPARM < 0.5, the number of iterations
required for convergence decreases dramatically as the ERRMAX criterion is
relaxed to values greater than 0.01 cm.
4.4 Transient Solution Strategy
The first step in establishing a transient solution strategy is to
estimate the time of leachate breakthrough using any of the available
analytical techniques. Three techniques reviewed in the Technical Resource
Document (KPA/530-SW-84-001) include the: (1) transit time equation, (2) the
11-29
-------
200
i
Oi
o
E
o
c
Q)
•4->
O
Q.
O
«•.
-100 -
D NODE 25
+ NODE 55
O NODE 85
10 13 16 19 22 25 28 31 34 37 40 43 46 49
-200 -
-300 -
-400
Number of Iterations
Figure 4-2. Steady state algorithm leads to numercial oscillation as shown at three selected nodes
from a clay/sand liner configuration (SRPARM = 1.0).
-------
NODE 85
200
100 -
£
o
c
-------
TABLE 4-1. EFFECT OF SRPARM ON THE NUMBER OF ITERATIONS
TO CONVERGENCE (90 cm LINER)a
HEAD (cm)
100.00
50.00
30.50
2.54
0.2
40
42
47
77
0.4
21
26
32
101b
SRPARM
0.5
17
26
35
101b
0.6
20
32
52
101b
0.8
101b
101b
101b
101b
a ERRMAX = 0.01 cm, MAXIT = 100
b MAXIT exceeded before convergence
TABLE 4-2. EFFECT OF SRPARM ON THE NUMBER OF ITERATIONS
TO CONVERGENCE (180 cm LINER)3
HEAD (cm)
100.00
50.00
30.50
2.54
0.2
44
43
50
83
0.4
23
34
45
101b
SRPARM
0.5
25
42
77
101b
0.6
43
101b
101b
101b
0.8
101b
101b
101b
101b
a ERRMAX =0.01 cm, MAXIT = 100
b MAXIT exceeded before convergence
11-32
-------
Green-Ampt wetting front model, and (3) the transient, linearized solution. A
first approximation of breakthrough serves as an aid in determining the total
duration (KNDTIM) of a transient simulation. Once an ENDTIM is established,
special output times can be selected to evaluate varying positions of the
wetting front with time. In many cases, however, steady state is achieved
before the specified ENDTIM. SOILINER terminates the time stepping procedure
and only tracks the particle based on the steady state PSI distribution.
Consequently, the last printed output time occurs at steady state, in which
case special output times greater than the time at which steady state is
achieved, and ENDTIM will not be printed. The remaining time stepping
parameters will be discussed after a brief description of the infiltration
event .
During the infiltration event, raatric potential, PSI, conductivity,
and moisture content, 6(40, are continuously changing until steady state is
achieved. This simultaneous change in soil properties as a function of PSI
effects the way in which PSI is calculated over time. To handle this
nonlinearity, SOILINER employs an iterative approach based on the finite
difference method. First, the new PSI distribution is predicted at a
calculated time level tj based on the previous PSI distribution and
associated soil properties. The soil properties are then updated at tj to
reflect the newly calculated PSI distribution. Based on the new soil
properties, SOILINER attempts to adjust each nodal value of pressure, PSI, by
a calculated change in pressure, DPSI.
Each time SOILINER updates the soil properties and readjusts PSI by the
value DSPI (see 3-15, Part I), an iteration is completed at the time level
tj. SOILINER continues to iterate until the maximum DPSI at all nodes
becomes less than the user-specified, error for convergence, ERRMAX. To
prevent uncontrolled iteration for a solution which is not converging at t^,
a maximum number of iterations per time level (MAXIT) is also specified. If
in attempting to reduce DPSI to a value smaller than ERRMAX, MAXIT is
exceeded, a forced exit from the time level occurs and a warning is issued.
In either case, whether ERRMAX is met or MAXIT is exceeded, the iterative
procedure terminates and a new time step is determined.
In order to step forward in time, SOILINER requires input data which is
used to calculate the size of each time step, DT. One of these input
parameters is DT itself, which describes the initial time step increment.
After the first time step using the initial, user-specified DT, SOILINER
calculates each subsequent time step as:
_ DT * CHPARM ...
DT
where CHPARM is a time step change parameter, input by the modeler, and CHMAX
represents the maximum change in matric potential at any node within the
solution domain, and is calculated between successive time steps in the
program. Therefore, after the first time step, SOILINER forces a small time
step if CHMAX is large, and a large time step if CHMAX is small. This
algorithm is a characteristic of the model designed to provide an internal
method of assuring model efficiency during the initial phase of infiltration.
As time progresses towards steady state, changes in the PSI distribution
11-33
-------
between successive time steps is generally small. As a result, calculated
time steps may get excessively large. To prevent this occurrence, SOILINER
incorporates another temporal input parameter (DTMAX), which regulates the
maximum size of a calculated time step. If the calculated time step is
greater than DTMAX, SOILINER forces DT to equal DTMAX. This parameter
prevents the time steps from getting so large that the predictions of PSI
become inaccurate, which may lead to difficulties in iteration convergence and
unusually high run times.
SOILINER continues to step through time until steady state or the final
simulation time step (ENDTIM) is reached. Each new time step is generated as
DPSI falls below ERRMAX, or MAXIT is exceeded. However, as another method of
avoiding uncontrolled CPU time, SOILINER requires that the user input a
maximum number of time steps per simulation (MAXNT). Model execution will be
terminated if MAXNT is exceeded.
The transient solution technique employed by SOILINER also incorporates a
temporal weighting parameter, ALPHA, which is used to weigh the approximations
of the state variable, PSI, at a point in time, tn+^, between the known
values of PSI at tn, and the unknown values of PSI at tn+1. ALPHA can
vary between 0.0 and 1.0. When ALPHA = 0.0, equation 3.19 (Part I) is solved
explicitly based on the values of PSI at tn, which are known. When
ALPHA = 1.0, equation 3.19 (Part I) is solved fully implicitly because the
values of PSI at tn+^ are unknown. Values of ALPHA between 0.0 and 1.0
weight the approximations accordingly. A value of 0.5 for ALPHA assumes that
the best approximation of PSI at tn is based equally on the known values
of PSI at tn and the unknown values of PSI at tn+1.
4.4.1 Method of Temporal Sensitivity Analysis—
Several transient simulations were conducted on an IBM 3033 mainframe in
which DT, DTMAX, CHPARM, and ALPHA were varied in order to evaluate SOILINER1S
sensitivity to these parameters. All remaining input data were held constant
including the initial and boundary conditions, soil properties, and grid
discretization. Each transient simulation was conducted such that the final
output time (ENDTIM) was specified to assure that the simulation reached
steady state conditions. Then a steady state solution was obtained for these
same conditions for comparison with the transient simulations to determine how
variations in DT, DTMAX, CHPARM, and ALPHA affected model stability, accuracy
and efficiency.
Each run's stability, or ability to converge, was evaluated based on a
review of the output. The accuracy of each simulation was determined by
quantitative comparisons of the matric potential at ENDTIM with the steady
state solution. Good comparisons were considered as those which showed no
greater than 1.0 cm of difference in matric potential between any two
corresponding nodes. Although some simulations were allowed extended CPU
times to observes the results, generally one minute was arbitrarily chosen as
the maximum amount of CPU time allocated for a given simulation. Results of
the sensitivity analyses substantiate that one minute is an economic and
realistic amount of time when compared with the values of ERRMAX and MAXIT
used for these analyses. The efficiency of each simulation was evaluated
11-34
-------
based primarily on the total number of iterations (NKITER) for a given
simulation. If job summaries are not available to obtain CPU times, the
variable NKITER can be utilized as an indicator of CPU time.
The scenario employed for all the sensitivity analyses is outlined below
and illustrated in Figure 2-1 of Part I. This scenario was used to minimize
the total CPU time required for the large number of simulations performed in
this sensitivity analysis.
1) Total solution domain, zw (Figure 2-1, Part I), of -361.00 cm,
corresponding to 76 nodes (or 75 elements).
2) Clay soil liner depth, d, of 61.00 cm, with a hydraulic conductivity
of 1.0x10"' cm/sec, porosity of 0.4950, and 56 nodes of variable
grid size.
3) Natural sand fill of depth 300.00 cm, with a hydraulic conductivity
of 9.444x10"^ cm/sec, porosity of 0.2870, and 20 nodes of constant
15.0 cm grid size.
4) A constant pressure boundary condition of 100.00 cm of impounded
liquid, hi, above the clay liner.
5) A constant pressure boundary condition of 0.0 cm for the water table
at the base of the sand fill (-361.00 cm).
6) Constant initial matric potential in the clay liner of -300.00 cm
and variable initial raatric potential in the sand fill to establish
static equilibrium.
7) Temporal parameters of DTMAX = 8.64xl05 sec, MAXNT = 2000,
MAXIT = 50, and ERRMAX =0.01 cm.
These values will provide some direction for users, however, the modeler
should recognize that each scenario will probably require at least some
variation from the above in order to meet specific requirements for
convergence, accuracy, and efficiency.
4.4.2 Choosing a Time Step—
The considerations in selecting time step sizes are analogous to those in
grid design. One efficient numerical scheme uses very small time steps at the
beginning of infiltration, when the matric potential gradients are highest,
and gradually increases the time step size as the moisture front advances into
the soil liner and disperses. Again, the initial time step size and the rate
of gradation are usually determined during initial simulations of the problem
at hand. Given that initial time steps may be as low as one second, or less,
a variable time step size is imperative to maintain cost effectiveness during
long-term simulation.
11-35
-------
Based on the algorithm used to calculate DT (equation 4-2), it is clear
that, after the first time step, DT can be regulated to an extent by CHPARM.
However, there are several additional input parameters of equal importance
which need to be considered in selecting a value of DT. These parameters
include: DTMAX, ERRMAX, MAXIT, and MAXNT. Various combinations of these
parameters will lead to different results for a given scenario, with varying
levels of accuracy and efficiency. In some instances, the model may not
converge to a solution. There is no unique combination of these parameters
for a given problem which will yield the best set of results. Rather, the
modeler must determine what level of accuracy is required and weigh this
against the cost (in CPU time) of obtaining those results.
Based on the scenario used for the sensitivity analyses, SOILINER
appeared very sensitive to values of CHPARM and DTMAX, and rather insensitive
to values of DT ranging from approximately 10 to 500. The most efficient and
accurate results obtained (for the given scenario) during the sensitivity
analyses occurred with 50 < CHPARM< 100 and DT within the range stated above.
However, it should be noted that high values for either DT or CHPARM should be
moderated by relatively low values for the other. As can be seen in
Figure 4-4, SOILINER was not sensitive to the initial DT when an appropriate
value of CHPARM was utilized. Changes in DT of more than an order of
magnitude do not appear to affect CPU time significantly (Table 4-3).
However, as DT was increased above 100 for the simulations shown in Table 4-3,
MAXIT warnings occurred and the solution became progressively less accurate.
In addition, values of DT that were excessively large when combined with high
values of CHPARM caused inaccurate predictions of PSI for the first time
step. The user is cautioned against setting excessively high DT values since
the results from each time step are used as the basis for future predictions.
If a large initial time step is used, it should be noted that model
convergence does not necessarily guarantee accuracy. Furthermore, convergence
may not be achieved, as depicted in Figure 4-5.
The user is also cautioned against small values of DT based on the
following logic. Results indicated that if the calculated time step was
forced to remain small using CHPARM, the number of iterations required for
convergence per time step was low. Since there was little change in the
predicted value of PSI and the soil properties at each new time step, DPSI
reached a value smaller than ERRMAX in relatively few iterations. However,
MAXNT was set at a large value because the simulation required many time steps
to reach even the first specified output time. If the specified simulation
time was very large relative to the calculated values of DT, the CPU time, and
therefore the cost of the simulation was high (see Table 4-4).
If the simulation progress through time normally with no execution
warnings (i.e., MAXIT exceeded), DT will generally increase as steady state is
approached. Successively larger DT values arise as the change in PSI between
time steps decreases (See equation 4-2). Even though changes in PSI can get
relatively small, it has been shown that DT can become excessive, leading to
model instability. To prevent unusually large, calculated DT values, the
variable DTMAX is specified. If DT exceeds DTMAX, DT is set equal to DTMAX
within SOILINER. However, it is still possible to specify a DTMAX that is too
large.
11-36
-------
E
o
15
*-M
c
"o
Q_
100
50 -
-5O -
-100 -
T - 1 5O -
-200 -
-25O H
-3OO
30
O
A
X
STEADY STATE
DT = 10 sees.
DT = 50 sees.
DT = 100 sees
DT = 500 sees
—r~
50
Depth (cm)
Figure 4-4. SOILINER sensitivity to initial DT values (pressure distribution at ENDTIM of transient
simulation where CHPARM = 75 and ALPHA = 1.0).
-------
TABLE 4-3. CPU TIME AS A FUNCTION OF DT FOR CONSTANT CHPARM=100*
DT CPU Time/Simulation (Sec)
10 24.09
100 23.73
1,000 23.74
10,000 20.29
86,400 18.65
* All remaining input parameters constant for each simulation
11-38
-------
15O
100
E
u
c
0)
•4-1
o
CL
D STEADY STATE
+ DT = 1000 sees.
CHPARM = 1000
O DT = 86,400 sees. (1 DAY)
CHPARM =500
Depth (cm)
-300
0
100
Figure 4-5. Divergence associated with large DT and CHPARM values.
-------
TABLE 4-4. SAMPLE CPU TIME NECESSARY FOR SIMULATIONS WITH
SMALL DT AND CHPARM VALUES*
ALPHA
0.75
1.00
1.00
DT
10
10
10
CHPARM
10
10
1.0
CPU TIME/ SIMULATION
62.28
65.43
53.38a
(SEC)
All remaining input parameters constant for each simulation
Forced exit at this time due to MAXNT exceeded. This run terminated at
simulation time 4.5 hours (ENDTIM was set to 1 year).
11-40
-------
Figure 4-6 depicts changes in PSI distributions overtime for a two-foot
liner. Note that special output time 5 (2 years) represents a time step in
which MAXIT was also exceeded (overall, MAXIT was exceeded 8 times). This
problem was easily rectified by reducing DTMAX from 100 to 10 days, whereby
all MAXIT warnings were eliminated. Although the number of time steps
increased from 164 to 232, NKITER decreased from 1860 to 1634, a reflection of
improved model efficiency.
4.4.3 Choosing a Temporal Weighting Parameter—
The explicit solution (ALPHA = 0.0) for PSI at a future time step is
based on the known values of PSI at the present time step. This solution
technique should only be used with small time steps, which reduces model
efficiency. Also, even with small time steps the explicit solution is
conditionally stable (i.e. convergence may not occur). Conversely, the
implicit solution (ALPHA = 1.0) is unconditionally stable (i.e. it always
converges), but since the solution for PSI at a future time step is based on
the unknown derivatives of PSI at that future time, the accuracy of the
results is still a function of the time step size, although not as sensitive
as the explicit solution.
In order to evaluate the effect various values of ALPHA had on the
solution to a given scenario, several simulations were conducted in which
ALPHA was set at 0.0, 0.25, 0.50, 0.75, and 1.00. For each value of ALPHA,
the temporal parameters were varied in an attempt to achieve model convergence
within 59 seconds of CPU time while noting accuracy of the results.
Several simulations were conducted for values of ALPHA = 0.75 and 1.00.
It was found that these simulations were unconditionally stable, however, the
user is cautioned to review the output carefully because accuracy was
dependent on the time stepping parameters. For both 0.75 and 1.00, each
simulation was more easily adjusted (than was possible for lower values of
ALPHA) so that no maximum iteration warnings were issued. These adjustments
were made primarily with CHPARM, DT, and DTMAX. As indicated by Table 4-5, it
appears that transient solutions are most accurately and efficiently simulated
when the fully implicit solution (ALPHA = 1.00) is employed. It should be
noted that for each run shown with MAXIT warnings in Table 4-5, the final
output time (ENDTIM) was in good comparison with the steady state solution.
Thus, MAXIT warnings do not necessarily indicate poor results, although an
ideal simulation would not give any. The modeler is advised that as ALPHA is
increased, the solution technique increasingly estimates the state variable,
PSI, at time level tn+^ based on the predicted derivatives of PSI at that
11-41
-------
£
o
c
0)
-t-1
o
Q.
100
50
0
-50 -
-100 -
Changes in Pressure Distribution
•c -150 -
-200 -
-250 -
-300
ERR = 159 cm
ERRMAX = 0.01 cm
0
Depth (cm)
80
Figure 4-6.
Results of MAXIT warning due to excessive DTMAX, where the magnitude of error between
successive iterations (with respect to ERRMAX) can be determined from ERR.
100
-------
TABLE 4-5. SOLUTION CHARACTERISTICS AS A FUNCTION OF TEMPORAL
PARAMETERS AND MAXIT
ALPHA
0.75
0.75
0.75
0.75
0.75
0.75
0.75
1.0
1.0
1.0
1.0
1.0
1.0
1.0
DT
10
10
10
10
1000
1000
1000
10
10
10
10
50
100
500
CHPARM
10
75
50
100
100
100
100
10
50
75
100
75
75
75
MAXIT
50
50
50
50
225
175
100
50
50
50
50
50
50
50
CPU TIME
62.28
24.90
26.82
23.56
26.39
25.33
23.91
65.40
30.78
28.14
28.28
28.00
28.09
27.39
NUMBER OF MAXIT
EXCEEDED WARNINGS
2
1
0
1
1
1
1
1
1
0
4
0
0
0
11-43
-------
time. However, for a scenario other than the one used in this sensitivity
analysis, it is possible that solution techniques which weigh the
approximations of PSI somewhat more on the known values of PSI at time level
tn (e.g., ALPHA < 0.75) may be more efficient and accurate.
Model convergence with accurate results could not be obtained for any
value of ALPHA < 0.50. Problems encountered in these simulations were
primarily register overflows or underflows, or insufficient CPU time
(i.e., greater than 59 seconds - Table 4-6). Those simulations that did not
converge showed highly inaccurate results. Some simulations revealed raatric
potentials within the liner that exceeded the boundary conditions. It appears
that both the Crank-Nicolson and explicit methods of solution to the
unsaturated, transient flow equation should not be used with SOILINER. This
limitation exists because explicit solutions are conditionally stable and need
to be used with small time steps. Furthermore, although Crank-Nicolson
solutions are unconditionally stable, their accuracy is also a function of the
time step size. Since SOILINER uses an internal method of determining each
time step, the user cannot assure that each calculated DT value will remain
within the range required for a particular solution technique. Explicit and
Crank-Nicolson solutions typically showed convergence and/or inaccuracy
problems and it is recommended that these solution strategies be avoided.
4.4.4 Setting a Tolerance for Iteration Convergence—
The CPU requirements should be borne in mind when the tolerance for
iteration convergence, ERRMAX, is set. Values of ERRMAX that are very low
will typically require high values for MAXIT and lengthy simulation times. In
some cases, if ERRMAX is too low (i.e., less than 0.01 cm), convergence may
never occur regardless of the specified MAXIT. Conversely, if ERRMAX is too
high, run time will be minimized (by relaxing the requirement on iteration
convergence) at the expense of model accuracy. Typical values for ERRMAX
range from 1.0 to 0.01 cm.
Normally, the modeler should decide what level of accuracy is acceptable
and set ERRMAX to this value. One should then attempt a simulation using a
conservative value of MAXIT and adjust DT and CHPARM at this MAXIT value while
noting CPU time and model accuracy. In this manner, the user can
systematically identify the optimum temporal input parameters. It is the
modeler's decision to increase MAXIT (and therefore the cost of a simulation)
based on the results of prior simulations.
4.4.5 Choosing the Maximum Number of Iterations Per Time Step—
The maximum number of iterations per time step, MAXIT, is used as one
means of limiting the amount of CPU time for a given scenario. Depending on
the scenario and input parameters, the number of iterations required per time
step for convergence may vary greatly. Under most conditions, MAXIT was found
not to have exceeded 50. However, the modeler should review the output from a
given scenario to determine the benefit of increasing MAXIT directly if
convergence appeared to be occurring when a forced exit from the time step
occurred. If the solution appeared to be diverging (as shown in Figure 4-6),
other parameters must be changed.
11-44
-------
TABLE 4-6. VALUES OF DT AND CHPARM PRODUCING REGISTER
UNDERFLOWS/OVERFLOWS OR ABENDS AS A FUNCTION
OF ALPHA
ALPHA
DT
CHPARM
COMMENT
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.25
0.25
0.25
0.25
0.50
0.50
0.50
0.50
0.5U
0.50
0.50
0.50
0.50
0.50
0.50
0.1
0.1
0.1
10.0
100.0
1,000.0
5,000.0
0.1
1.0
10.0
10.0
0.1
0.1
0.1
1.0
1.0
10.0
10.0
10.0
100.0
1,000.0
10,000.0
0.1
10.0
100.0
10.0
100.0
1,000.0
5,000.0
0.1
1.0
10.0
1,000.0
1.0
50.0
1,000.0
1.0
50.0
1.0
10.0
1,000.0
50.0
50.0
50.0
A!
u/o2
u/o
u/o
u/o
u/o
u/o
A 2 min.
A 2 min.
A 2 rain.
U/O 2 rain.
U/O
U/O
A
A
U/O
U/O
U/O
A
U/O
U/O
U/O
CPU allocated
CPU allocated
CPU allocated
CPU allocated
Abend (i.e., Program terminated due to run-time errors or insufficient
CPU time.)
Underflow or Overflow (i.e., Numbers were too small or large,
respectively, for the computer to handle.)
All remaining input parameters constant for each simulation
All simulations allocated 50 CPU seconds unless otherwise noted.
11-45
-------
In some cases, the impact from exceeding MAXIT may not be as obvious as
that shown in Figure 4-6. For the last time step shown in Figure 4-6 the
calculated ERR at MAXIT was 159.5 cm, far greater than the specified ERRMAX of
0.01 cm. In other cases, it is possible that the criterion established for
convergence was barely missed (i.e., ERR = 0.02 cm with an ERRMAX set at
0.01 cm) when the MAXIT warning was issued. For this situation, the
calculated PSI distribution associated with the MAXIT warnings may be
reasonable. Since MAXIT warnings are frequently issued for a time step that
does not coincide with a special output time, it is difficult to visualize the
magnitude of the ERR associated with the warning. When MAXIT is exceeded,
erroneous values for PSI and the associated soil properties are used as the
basis for subsequent predictions through the forcing vector F, and stiffness
matrix STIFF (see Part I and Appendix F). Consequently, an extremely poor
prediction at a given time step (Figure 4-6) may be propagated through time
even though subsequent time steps may converge with no MAXIT warnings. As an
aid in determining the impact of MAXIT warnings on the overall simulation,
SOILINER provides the calculated ERR each time a MAXIT warning is issued.
Ideally, the parameters will be adjusted accordingly such that no warnings are
issued for a desired simulation.
4.4.6 Sensitivity of the Particle Tracking Algorithm—
Since the particle tracking algorithm is based on the predicted PSI
distribution at each time step, it is very sensitive to prediction accuracy.
Figure 4-7 demonstrates, by comparison with Figure 2-5, the impact of model
instability on particle movement. The results depicted in both graphics were
obtained from the same simulation, with one exception. The variable, DTMAX,
was changed from 10 days (Figure 2-5) to 100 days as shown by the increased
spacing between data points towards the latter half of the simulation time in
Figure 4-7. The results in this figure were derived from the predicted PSI
distributions of Figure 4-6. For the time steps numbered in Figure 4-7, MAXIT
was exceeded and erroneous PSI results were obtained and in some instances,
caused the particle to move upward into the liner. By readjusting parameters
(i.e., DTMAX reduced to 10 days) it was possible to eliminate all warning
messages resulting in an accurate simulation of particle movement (Figure 2-5),
11-46
-------
E
u
Q.
O
Q
0
-10 -
-20 -
-30 -
-40 -
-50 -
-60
0
Particle Location (61cm Clay Liner)
0.8
1.2 1.6
Time (yrs.)
2.4 2.8
Figure 4-7.
Instability of particle movement due to inaccurate
PSI distributions associated with divergence at the
time steps shown above.
-------
REFERENCES
Anderson, D.C., and S.G. Jones, Clay barrier-leachate interaction, presented
at National Conference on Management of Uncontrolled Hazardous Waste
Sites, Washington, DC, October 31-November 2, 1983.
Bear, R.J., Hydraulics of Groundwater, McGraw-Hill, New York, 1979.
Brooks, R.H., and A.T. Corey, Hydraulic properties of porous media,
Hydrology Paper No. 3, Colorado State University, Fort Collins,
March 1964.
Brutsaert, W.F., A function iteration technique for solving the Richards
equation applied to two-dimensional infiltration problems, Water
Resources Research 6(7), 1583-1596, 1971.
Campbell, G.S., A simple method for determining unsaturated conductivity from
moisture retention data, Soil Sci., 117, 311-314, 1974.
Clapp, R.B., and G.M. Hornberger, Empirical equations for some soil hydraulic
properties, Water resources Research, 14, 601-604, August 1978.
Cooley, R.L., A finite difference method for unsteady flow in variably
saturated porous media: application to a single pumping well, Water
Resources Research, 7(6), 1607-1625, December 1971.
Cooper, H.H., Jr., J.D. Bredehoeft, and I.S. Papadopoulos, Response of a
finite-diameter well to an instantaneous charge of water, Water Resources
Research, 3, 263-269, 1967.
Dane, J.H., and P.J. Wierenga, Effect of hysteresis on the prediction of
infiltration, redistribution, and drainage of water in a layered soil,
Journal of Hydrology, 25, 229-242, 1975.
Elzeftawy, A., and K. Cartwright, Evaluating the saturated and unsaturated
hydraulic conductivity of soils, in T.F. Zimmie and C.O. Riggs, eds.,
Permeability and Groundwater Contaminant Transport, American Society for
Testing and Materials, 168-181, 1981.
Elzeftawy, A., and B.J. Dempsey, Unsaturated transient and steady state flow
of moisture in subgrade soil, Transportation Research Rec. 612, 56-61,
1976.
11-48
-------
Freeze, R.A., The mechanism of natural ground-water recharge and discharge:
1. One-dimensional, vertical, unsteady, unsaturated flow above a
recharging or discharging ground-water flow system, Water Resources
Research, 5(1), 153-171, 1969.
Freeze, R.A., and A. Cherry, Groundwater, Prentice-Hall,
Englewood Cliffs, NJ, 1979.
Gardner, W.R. , Some steady-state solutions of the unsaturated moisture flow
equation with application to evaporation from a water table, Soil
Science, 8, 228-232, 1958.
Gillhara, R.W, A. Klute, and D.F. Heerman, Hydraulic properties of a porous
medium measurement and empirical representation, Soil Science Society of
America Journal, 40, 203-207, March-April 1976.
Goode, D.J., and P.A. Smith, Procedures for modeling flow through clay liners
to determine required liner thickness, draft technical resource document
for public comment, EPA/530-SW-84-001, 1984.
Green, D.W. , H. Dabiri, and C.F. Weinaug, Numerical modeling of unsaturated
groundwater flow and comparison of the model to a field experiment, Water
Resources Research, 6, 862-874, June 1970.
Haverkamp, R., M. Vauclin, J. Toma, P.J. Wierenga, and G. Vachaud, A comparison
of numerical simulation models for one-dimensional infiltration, Soil
Science Society of America Journal, 41, 285-294, 1977.
Hillel, D.I., Soil and Water - Physical Principles and Processes, Academic
Press, New York, 1971.
Hillel, D.I., Fundamentals of Soil Physics, Academic Press, N.Y., 1980.
Hvorslev, M.J., Time lag and soil permeability in groundwater observations.
U.S. Army Corps of Engineers Waterways Exp. Sta. Bull. 36, Vicksburg,
Miss, 1951.
Jackson, R.D., R.J. Reginato, and C.H.M. Van Bavel, Comparison of measured and
calculated hydraulic conductivities of unsaturated soils, Water Resources
Research, 1(3), 375-380, 1965.
Johnson, R.A., and E.S. Wood, Unsaturated flow through clay liners, letter
report in fulfillment of EPA Contract No. 68-01-6871, GCA/Technology
Division, Bedford, MA, 1984.
King, L.G., Description of soil characteristics for partially saturated flow,
Soil Science Society of America Proceedings, 29(4), 359-362, 1965.
Kunze, R.J., and D.R. Nielsen, Finite-difference solutions of the infiltration
equation, Soil Science, 134(2), 81-88, August 1982.
11-49
-------
Li, E.A., V.D. Shanholtz, and E.W. Carson, Estimating saturated hydraulic
conductivity and capillary potential at the wetting front, Dept. of Agr.
Eng., VA Polytech. Inst. and State Univ., Blacksburg, 1976.
McQueen, I.S., and R.F. Miller, Approximating soil moisture characteristics
from limited data: empirical evidence and tentative model, Water
Resources Research, 10(3), 521-527, June 1974.
Milly, P.C.D., and P.S. Eagleson, The coupled transport of water and heat in a
vertical soil column under atmospheric excitation, R.M. Parsons
Laboratory for Water Resources and Hydrodynamics, M.I.T., Technical
Report No. 258, July 1980.
Morel-Seytoux, H.J., Two-phase flows in porous media, 119-120, in V.T. Chin,
ed., Advances in Hydroscience, Academic Press, N.Y., 1973.
Mualem, Y., A new model for predicting the hydraulic conductivity of
unsaturated porous media, Water Resources Research, 12(3), 513-522,
June 1976.
Mualem, Y., Hydraulic conductivity of unsaturated porous media: generalized
microscopic approach, Water Resources Research 14(2), 325-334, April 1978.
Philip, J.R., The theory of infiltration: 6. Effect of water depth over soil,
Soil Science, 85(5), 278-286, 1958.
Philip, J.R., Theory of infiltration, 215-297, in V.T. Chow, ed., Advances in
Hydroscience, Vol 5, Academic Press, N.Y., 1969.
Pinder, G.F., and W.G. Gray, Finite Element Simulation in Surface and Sub-
surface Hydrology, Academic Press, N.Y., 1977.
Papadopoulos, I.S., J.D. Bredehoeft, and H.H. Cooper, On the analysis of slug
test data, Water Resources Research, 9, 1087-1089, 1973.
Ragab, R., J. Feyen, and D. Hillel, Comparison of experimental and simulated
infiltration profiles in sand, Soil Science, 133(1), 61-64, January 1982.
Reeder, J.W., D.L. Freyberg, J.B. Franzini, and I. Remson, Infiltration under
rapidly varying surface water depths, Water Resources Research, 16(1),
97-114, February 1980.
Rogowski, A.S., Watershed physics: model of the soil moisture characteristic,
Water Resources Research, 7(6), 1575-1582, December 1971.
11-50
-------
APPENDIX A
LIST OF SYMBOLS USED IN TEXT
A-l
-------
Symbol
9010
e
r
15
k
K
K
FORTRAN
Variable
ALPHA
RMOIST
PSI
PSICRT
DP SI
C, COLD
H
RELK
SATK
RKL
POR
Definition
Temporal weighting parameter
Specific weight of water
Moisture content (equals n at saturation
point)
Average element moisture content
Non-reducible moisture content
Saturated water content
Moisture content at 1.5x10^ cm pressure
Total head (
-------
FORTRAN
Symbol Variable Definition
q Darcy flux
qe FLUX Darcy element flux
Se Effective saturation
t TIME Time
At DT Change in time (length of time step)
V(40 VNEW, VOLD Velocity vector (RHS of equation 3-9,
Part I)
ve VELO Element interstitial velocity
W Soil wetness
z Z Depth below datum at liner surface
Az DZ Distance between two adjacent nodes
zw Depth below datum of the water table
A-3
-------
APPENDIX B
COMPUTER PROGRAM FOR GARDNER'S ANALYTICAL SOLUTION
B-l
-------
09/26/83 08:40:42 GCA6D.GARDNER.FORT.DATA
K =
00000310 C«*«
00000020 C
OOOOOr-30 C
000000*0 C
00000050 C
00000060 C
00000070 C
00000080 C
OQOOOQ9Q C
00000100 C
00000110 C
00000120 C
00000130 C
00000140 C
00000150 C
00000160 C
00000170 C
00000180 C
00000190 C*••••••••••»•»»•
00000200 C
00000210
00000220
00000230
00000240 C
00000250
00000260
00000270 C
00000280
00000290
00000300
00000310
00000220
00000330
00000340
00000350
00000360
00000370
00000380 C
00000390 C
00000400
00000410
00000420 C
00000430
00000440 C
00000450
00000460
00000470
C0000480 C
00000490
00000500
00000510
00000520
00000530
00000540
00000550 C
00000560
00000570
00000580
00000590 C
00000600
GARDNER.FORT STEADY STATE EVAPORATION IN UNSATURATEu SOIL
SEE y.R. GARDNER* SOME STEADY-STATE SOLUTIONS OF THE UNSATURATED
MOISTURE FLOW EQUATION ylTH APPLICATION TO EVAPORATION FROM
A WATER TABLE* SOIL SCIENCE* 85* 228-232* 1958.
STEADY STATE SOLUTION FOR SOIL WITH HYDRAULIC CONDUCTIVITY
FUNCTION AS FOLLOWS:
S = - (PSI) SUCTION
S * BETA
OAN GOOOE JUNE 1983
IRD=5
IPRT=6
IRSLT=9
READ SOIL FUNCTION PARAMETERS
READURD,*) A.B«WT.FLUX
WRITE(IPRT*2001> A,B*nT*FLUX
COMPUTE SOLUTION PARAMETERS
SQR2=1.4142136
ALPHA=FLUX/A
BETA=ALPHA*b « 1.0
R0=
T£RM2=1./(2.*R03»SQR2>
URITE«IPRT,2004) ALPHA.BETA,RO,R02.R03.AINV«TER«1.TERH2
READ STEPPING FOR SOLUTION
READtIRD**) PF.DPF.PFENO*MAX
URITEUPRT.2002) PF,DPF*PFENO*MAX
PF=PF-DPF
DO 10 1=1,MAX
PF=PF»OPF
IF(PF.GT.PFENO> GOTO 999
S=10.0**PF
PSIs-S
S2=S*S
TERM3:RO*S*SOR2
ARG1=CS2»TERM3»R02)/
ARG2=TERH3/(R02-S2)
PARH=ATAN«ARG2>
IF (ARG2.LT.O) PARM=PARM « 3.1415927
Z=(ALOG(ARG1)*TERM1*PARM*TERM2)*AINV-UT
URITE(IPRT«2003I I,Z,PSI»PF*S.S2»TERM3.ARG1.ARG2
B-2
-------
00000610
OOOCOt20
00000630
00000640
00000650
00000660
00000*70
00000660
00000690
00000700
30000710
30000720
OC000730
000007*0
OOOOC750
00000160
00000770
00000760
00000790
00000800
OOOOOblO
00000620
00000030
00000640
00000650
00000860
00000870
30000683
10
URITEIIRSLT,20031
CONTINUE
I.Z.PSI.PF
•999 STOP
2C01 FORMATC//' GARDNER OUTPUT'//
••.SOIL FUNCTION PARAMETERS' //
•* A'.lSdH.) t'A = »,lPE10.2/
•• 6',15<1H.) »'B=',1PE10.2/
•• CE'TH TO HATER TABLE* .5J1H. >.'UT=',OPF10.2/
*• FLJX UPWARD'»7<1H.).'FLUX=»,1PE10.2>
2002 FORMAT!//' SOLUTION STEPPING PARAMETERS'//
*• BEGINNING SUCTION' « 10 < 1H. ). 'PF=' . OPF10 . 21
• • PF INCREMENT'. 10UH.),»DPF = ',OPF10.3/
*' FINAL PF'.12<1H.),'PFENO='.OPF10.2/
•» MAXIMUM NUMBER OF P 3INTS' .5 C1H. ) . 'MAX=' . 110///
*• I Z PSI PF S
•3 ARG1 ARG2',/)
2003 FORMAT(I10«GF>F10.3.1PE10.2,OPF10.3.6<1P£10.2»
2004 FOKMAT
-------
APPENDIX C
PARTIAL LIST OF FORTRAN VARIABLES
IN SOILINER
C-l
-------
ALPHAS - Temporary storage of the weighting parameter, ALPHA.
AMI - 1.0 - ALPHA.
AVEL - Average particle velocity of three successive time steps.
C - Array containing the average node values of the specific
moisture capacity, based on data from two adjacent elements
(see CL).
CDIST - Distance traveled by the particle from its last position to the
point of breakthrough.
CHANGE - Absolute value of the change in pressure at a given node
between successive time steps.
CL - Each interior, mesh centered node is broken into two nodes, one
associated with each adjacent element; the CL array contains
values of the specific moisture capacity (i.e., the derivative
of moisture with respect to pressure) for this expanded set of
nodes.
DELZ - Absolute value of the distance between two adjacent nodes.
DPSI - Array holding nodal values of the change in pressure between
two successive iterations.
DSTOR - Total moisture storage (cm) during time DT for all elements in
the flow domain.
DTV - Time required for the particle to travel from its previous
position to the point of breakthrough.
DZ - Distance (cm) between two adjacent nodes.
DZN - Distance between midpoint of two adjacent elements.
EFLUX - Mass balance absolute error (cm/sec) considering fluxes and
storage rate at a given time step.
EREL - Relative mass balance storage error.
ERR - The maximum, absolute value at a node for the change in
pressure between successive iterations.
EVOL - Mass balance absolute error (cm) considering input, output, and
storage per unit surface area during time step DT.
ETOT - Mass balance absolute error (cm) considering cutnraulative input,
ouput, and storage per unit surface area.
C-2
-------
F - Forcing vector of the matrix equation solved by the Thomas
algorithm (contains all knovms on the right hand side of
equation 3-19, Part I).
FLUXl - Flux into top element (cm/sec) of the flow domain at the
current time step.
FLUX10 - Fluxl from previous time step.
FLUX2 - Flux leaving the flow domain (cm/sec) at the current time step.
FLUX20 - Flux2 from previous time step.
IARG - Variable passed as an integer argument to subroutines.
IOUT - Integer flag indicating special output time to be printed.
HER - Number of k-iterations for a given time step.
IWARN - Total number of times maximum iteration warning is given.
JSOIL - Variable containing soil code number; used to determine
appropriate characteristic curve for a given element.
KODE - Integer flag used to indicate that information associated with
a special output time is to be printed.
KODEBT - Integer flag indicating time step when breakthrough occurred.
KODESS - Integer flag indicating time step when steady state achieved.
NELEMV - Element of particle location.
NKITER - Total number of k-iterations (number of times matrix equation
solved).
NMAX - Maximum number of nodes (set at 200).
NODEV1 - Upper node of NELEMV.
NODEV2 - Lower node of NELEMV.
NOUTl - Index used to determine which special output time is to be
printed next.
NPM1 - Number of nodes minus one.
NPM2 - Number of nodes minus two.
NT - Total number of time steps.
NUMEL - Total number of elements.
NUMEL.2 - Two times the number of time elements.
C-3
-------
NX - Number of data points for the velocity integration scheme using
Simpson's Rule.
OLDC - Array containing nodal values of the specific moisture capacity
from the previous time step.
OLDML - Array holding nodal values of moisture content from the
previous time step, for the expanded node set (see CL).
OLDTD - Total distance traveled by the particle at the previous time
step.
PSILN - Log Base 10 of the negative matric potential for a particular
node (PSI must be less than 0).
PSINEG - Negative value of PSI.
PSIOLD - Array holding pressure distribution from previous time step.
RATES! - Rate of liquid storage during infiltration over the period UT.
RELK - Fraction used to multiply the saturated conductivity value for
a given node to determine the relative conductivity of the node
when it is unsaturated.
RK - The calculated, unsaturated conductivity at a given node.
RKL - Array holding nodal values of conductivity for the expanded
node set (see CL).
RMOIST - Moisture content at a given node.
RMSTL - Array holding nodal values of moisture content for the expanded
set of nodes (see CL).
SRPSI - Array holding the pressure distribution as calculated using the
successive relaxation algorithm.
SSPSI - Array holding the final steady state pressure distribution.
STARK - Element hydraulic conductivity.
STOR - Cumulative moisture storage per unit surface are (cm) for all
elements in the flow domain up to the current time.
SYEAR - Seconds per year.
TD1ST - Total distance traveled by the particle.
C-4
-------
TIME1 - Temporary storage for the variable TIME.
TIMEBT - Time at the point of particle breakthrough (cm).
TOTV1 - Cumulative volume flux per unit area (cm) through the top
element of the flow domain.
TOTV2 - Cumulative volume flux per unit area (cm) through the bottom of
the flow domain.
TOUT1 - Variable holding the current special output time (after each
special output time, TOUT1 is assigned the next special output
time).
VEL - Calculated interstitial velocity (subroutine OUTPUT).
VELO - Array holding interstitial velocity data necessary for
integration of the velocity function.
VNEW - Array holding current values of the velocity vector
(equation 3-19, Part I) which is updated each interation.
VOLD - Array holding valves of the velocity vector (equation 3-19,
Part I) from the previous time step.
VOL1 - Volume flux per unit surface area (cm) through the top element
over the time period DT.
VOL2 - Volume flux per time unit surface area (cm) through the bottom
element over the time period DT.
VTIME - Array holding times required for integration of the velocity
function.
WORK - Working array containing coefficients of the stiffness matrix
used in the Thomas algorithm.
YDTV - DTV (years).
YTIME - Time in years.
ZDIST - Particle distance traveled over the time period of integration.
ZZ - Midpoint of an element.
C-5
-------
APPENDIX D
EXAMPLE OF INPUT DATA SETS
D-l
-------
This Appendix provides annotated copies of the four data sets required
as input to the SOILiNER model. These data pertain to a double liner system as
follows:
Layer 1: 61 cm (2") clay; K = 1.00 x 10~ cm/sec, n = 0.495
Layer 2: 30 cm (I1) sand; K = 1.76 x 10~ cm/sec, n = 0.395
Layer 3: 91 cm (3') clay; K = 1.00 x 10~7 cm/sec, n = 0.495
_2
Layer 4: 300 cm to the water table; K = 1.76 x 10 cm/sec, n = 0.395
0-2
-------
TEST: DOUBLE LINER; 61/30/91/300CM CLAY/SAND/CLAY/SANDj TRANSIENT
T F 4.0 100. 6.05E+5 3000 1.0 .10 50 80.
4
0.5 1.0 2.0 3.0
152 0.40 132
100.0 0.0
Control File
RECORD TYPE 1: FORMAT(80A)
TITLE - character string identi-fying run
RECORD TYPE 2: UNFORMATTED
LPRINT - boolean flag for output format (F « shortj T « long)
LSTEDY - boolean flag (F - transient solution; T = steady state)
ENDTIM - estimated solution time (yrs)
DT - initial time step size (sec)
DTMAX - maximum allowable time step size (sec)
MAXNT - maximum number of time steps
ALPHA - temporal weighting parameter (0.0 s ALPHA 5 1.0)
ERRMAX - convergence criteria (cm)
MAXIT - maximum number of iterations at any given time step
CHPARM - change parameter applied to step size between steps
RECORD TYPE 3: UNFORMATTED
NOUT - number of special output times (4 maximum)
RECORD TYPE 4: UNFORMATTED
** NOTE: If NOUT < = 0, this record must be excluded.
TOUT(1:NOUT) - special output times (yrs)
RECORD TYPE 5: UNFORMATTED
NUMNP - number of nodes (200 maximum)
SRPARM - successive relaxation parameter (0.0 s SRPARM s 1.0)
NCLAYN - number of liner nodes
RECORD TYPE 6: UNFORMATTED
H - head in impoundment (cm)
PSIDOT - pressure at base of flow domain (cm)
D-3
-------
10 0.10 Brid File
IB 0.50
10 1.00 RECORD TYPE 1: UNFORMATTED
10 1.50
15 2.00 MUM - number of consecutive nodes
6 5.00 (downwards) with separation DELTA
10 0.10 DELTA - separation between adjacent nodes
10 0.50
10 1.00
10 1.50
30 2.00
20 15.00
D-4
-------
131 -150. Initial Conditions File
1 -3(90.
1 -285. RECORD TYPE li UNFORMATTED
1 -270.
1 -255. NUM — number o-f consecutive nodes with
1 -240. pressure PS I
1 -225. PSI - initial pressure (cm)
1 -210.
1 -195.
1 -180.
1 -165.
1 -150.
1 -135.
1 -120.
1 -105.
1 -90.
1 -75.
1 -60.
1 -45.
1 -30.
1 -15.
1 0.
D-5
-------
55
6
713
20
12 1.000E-07
1 1.760E-02
12 1.000E-07
1 1.760E-02
0.4950
0.3950
0.4950
0.3950
-1.000
-16.96
-1.000
-16.96
Soil Properties File
RECORD TYPE 1: UNFORMATTED
NUM - number o-f nodes in soil layer
ISOIL - soil code number
SATK - saturated conductivity for soil ISOIL
POR - porosity -for soil ISOIL
PSICRT - critical pressure -for soil ISOIL
D-6
-------
APPENDIX E
EXAMPLE OF OUTPUT DATA
E-l
-------
This Appendix provides example output for the input data shown in
Appendix D. The first half of the output is primarily an echo of the input
(LPRINT = .TRUE.). Echoed data is followed by a brief summary of the results
from the steady state algorithm, including the number of iterations and the
error at convergence. If convergence had not been achieved, a warning would
have been issued.
A total of six output times is provided - initial conditions, the four
specified output times, and the steady state solution. After each output time
a mass balance summary is printed. Finally, both the times to steady state
and breakthrough are given when they occur. For this example steady state was
achieved at year 3.8 and breakthrough occured at year 12.0.
Only the General Output file is included in this Appendix. The four
remaining output files contain data in a format suitable for graphical
analysis. These files and their output data are summarized below:
Flux Output (FLX.OUT)
L - element number
ZZ - midpoint of the element (cm)
FLUX - element flux (cm/sec)
VEL - element interstitial velocity (cm/sec)
Pressure Output (PSI.PRN)
I - node number
Z - node depth (cm)
PSI - calculated pressure (cm)
PF - log base 10 of the nodal pressure
Moisture Output (MST.PRN)
I - node number
Z - node depth (cm)
RMOIST - moisture value at the node
Particle Depth Output (PDT.OUT)
NT - time step number
TIME - time (yrs)
TDIST - total distance travel by the particle since initiation of the
simulation (cm)
E-2
-------
SOILINER OUTPUT
TEST: DOUBLE LINER; 61730/91/300CM CLAY/SAND/CLAY/SAND; TRANSIENT
TEMPORAL DISCRETIZATION PARAMETERS
STEADY STATE PARM LSTEDY= F
IF LSTEDY EQ T, COMPUTE STEADY STATE ONLY
OTHERWISE, COMPUTE TRANSIENT SOLUTION
SIMULATION TIME ENDTIM(YRS)« 4.00
TIME STEP DT= . 100E+03
MAXIMUM ALLOWABLE TIME STEP DTMAX= .6050E+06
MAXIMUM NUMBER OF TIME STEPS MAXNT= 3000
TEMPORAL WEIGHTING PARAMETER ALPHA= 1.00
MAXIMUM ERROR FOR CONVERGENCE ERRMAX= .1000
MAXIMUM ITERATIONS PER TIME STEP MAXIT= 50
TIME STEP CHANGE PARAMETER CHPARM= 80.0000
NUMBER OF SPECIAL OUTPUT TIMES NOUT- 4
SPECIAL OUTPUT TIMES(YRS)
.50 1.00 2.00 3.00
SPATIAL DISCRETIZATION PARAMATERS
NUMBER OF NODE POINTS NUMNP= 152
NUMBER OF ELEMENTS (COMPUTED) NUMEL= 151
SUCCESSIVE RELAXATION PARAMETER SRPARM- .40
NUMBER OF LINER NODES NCLAYN* 132
(USED FOR BREAKTHROUGH DETERMINATION)
GRID DATA
NODE Z DZN
ELEMENT DZ
E-3
-------
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
1 .000
-.100
2 -.100
-.100
3 -.200
-.100
4 -.300
-.100
5 -.400
-.100
6 -.500
-.100
7 -.600
-.100
8 -.700
-.100
9 -.800
-.100
10 -.900
-.100
11 -1.000
-.500
12 -1.500
-.500
13 -2.000
-.500
14 -2.500
-.500
15 -3.000
-.500
16 -3.500
-.500
17 -4.000
-.500
18 -4.500
-.500
19 -5.000
-.500
20 -5.500
-.500
21 -6.000
-1.000
22 -7.000
-1.000
23 -8.000
-1.000
24 -9.000
-1.000
25 -10.000
-1.000
26 -11.000
-1.000
-.100
-.100
-.1B0
-.100
-.100
-.100
-. 100
-.100
-.100
-.300
-.500
-.500
-.500
-.500
-.500
-.500
-.500
-.500
-.500
-.750
-1.000
-1.000
-1.000
-1.000
-1.000
E-4
-------
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
27 -12.000
-1.000
28 -13.000
-1.000
29 -14.000
-1.000
30 -15.000
-1.000
31 -16.000
-1.500
32 -17.500
-1.500
33 -19.000
-1.500
34 -20.500
-1.500
35 -22.000
-1.500
36 -23.500
-1.500
37 -25.000
-1.500
38 -26.500
-1.500
39 -28.000
-1.500
40 -29.500
-1.500
41 -31.000
-2.000
42 -33.000
-2.000
43 -35.000
-2.000
44 -37.000
-2.000
45 -39.000
-2.000
46 -41.000
-2.000
47 -43.000
-2.000
48 -45.000
-2.000
49 -47.000
-2.000
50 -49.000
-2.000
51 -51.000
-2.000
52 -53.000
-2.000
-1.000
-1.000
-1.000
-1.000
-1.250
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.750
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
E-5
-------
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
53 -55.000
-2.000
54 -57.000
-2.000
55 -59.000
-2.000
56 -61.000
-5.000
57 -66.000
-5.000
58 -71.000
-5.000
59 -76.000
-5.000
60 -81.000
-5.000
61 -86.000
-5.000
62 -91.000
-. 100
63 -91.100
-.100
64 -91.200
-.100
65 -91.300
-. 100
66 -91.400
-.100
67 -91.500
-.100
68 -91.600
-. 100
69 -91.700
-.100
70 -91.800
-.100
71 -91.900
-. 100
72 -92.000
-.500
73 -92.500
-.500
74 -93.000
-.500
75 -93.500
-.500
76 -94.000
-.500
77 -94.500
-.500
78 -95.000
-.500
-2.000
-2.000
-2.000
-3.500
-5.000
-5.000
-5.000
-5.000
-5.000
-2.550
-. 100
-. 100
-.100
-.100
-.100
-.100
-.100
-.100
-. 100
-.300
-.500
-.500
-.500
-.500
-.500
-.500
E-6
-------
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
79 -95.500
-.500
80 -96.000
-.500
81 -96.500
-.500
82 -97.000
-1.000
83 -98.000
-1.000
84 -99.000
-1.000
85 -100.000
-1.000
86 -101.000
-1.000
87 -102.000
-1.000
88 -103.000
-1.000
89 -104.000
-1.000
90 -105.000
-1.000
91 -106.000
-1.000
92 -107.000
-1.500
93 -108.500
-1.500
94 -110.000
-1.500
95 -111.500
-1.500
96 -113.000
-1.500
97 -114.500
-1.500
98 -116.000
-1.500
99 -117.500
-1.500
100 -119.000
-1.500
101 -120.500
-1.500
102 -122.000
-2.000
103 -124.000
-2.000
104 -126.000
-2.000
-.500
-.500
-.500
-.750
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.250
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.500
-1.750
-2.000
-2.000
E-7
-------
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
105 -128.000
-2.000
106 -130.000
-2.000
107 -132.000
-2.000
108 -134.000
-2.000
109 -136.000
-2.000
110 -138.000
-2.000
111 -140.000
-2.000
112 -142.000
-2.000
113 -144.000
-2.000
114 -146.000
-2.000
115 -148.000
-2.000
116 -150.000
-2.000
117 -152.000
-2.000
118 -154.000
-2.000
119 -156.000
-2.000
120 -158.000
-2.000
121 -160.000
-2.000
122 -162.000
-2.000
123 -164.000
-2.000
124 -166.000
-2.000
125 -168.000
-2.000
126 -170.000
-2.000
127 -172.000
-2.000
128 -174.000
-2.000
129 -176.000
-2.000
130 -178.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
-2.000
E-8
-------
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
)IL
131 -180.000
-2.000
132 -182.000
-15.000
133 -197.000
-15.000
134 -212.000
-15.000
135 -227.000
-15.000
136 -242.000
_ i >5 nnm
i^J m KIKIKI
137 -257.000
-15.000
138 -272.000
-15.000
139 -287.000
-15.000
140 -302.000
-15.000
141 -317.000
-15.000
142 -332.000
-15.000
143 -347.000
-15.000
144 -362.000
-15.000
145 -377.000
-15.000
146 -392.000
-15.000
147 -407.000
-15.000
148 -422.000
-15.000
149 -437.000
-15.000
150 -452.000
-15.000
151 -467.000
-15.000
152 -482.000
PROPERTIES
ELEMENT I SOIL
1 12
2 12
3 12
-2.000
-8.500
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
-15.000
SATK
1 . 000E-07
1 . 000E-07
1 . 000E-07
POR
.4950
.4950
.4950
PSICRT
-1.000
-1.000
-1.000
E-9
-------
4
5
6
7
8
9
10
11
12
13
14
15
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
48
49
50
51
52
53
54
55
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
12
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
1 . 000E-07
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
E-10
-------
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
.3950
.3950
.3950
.3950
.3950
.3950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
E-ll
-------
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
12 1.000E-07
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1 . 760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
1 1.760E-02
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
.3950
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
-16.960
INITIAL PSI
NODE
1
2
PSI
-150.00
-150.00
E-12
-------
3
4
5
6
7
8
9
10
11
12
13
14
15
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
48
49
50
51
52
53
54
•150.00
•150.00
-150.00
-150.00
-150.00
•150.00
-150.00
•150.00
-150.00
•150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
E-13
-------
55
56
57
SB
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
-150.00
•150.00
-150.00
•150.00
-150.00
•150.00
-150.00
-150.00
-150.00
•150.00
-150.00
•150.00
-150.00
-150.00
-150.00
•150.00
-150.00
•150.00
-150.00
•150.00
-150.00
-150.00
-150.00
•150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
•150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
-150.00
E-14
-------
107 -150.00
108 -150.00
109 -150.00
110 -150.00
111 -150.00
112 -150.00
113 -150.00
114 -150.00
115 -150.00
116 -150.00
117 -150.00
118 -150.00
119 -150.00
120 -150.00
121 -150.00
122 -150.00
123 -150.00
124 -150.00
125 -150.00
126 -150.00
127 -150.00
128 -150.00
129 -150.00
130 -150.00
131 -150.00
132 -300.00
133 -285.00
134 -270.00
135 -255.00
136 -240.00
137 -225.00
138 -210.00
139 -195.00
140 -180.00
141 -165.00
142 -150.00
143 -135.00
144 -120.00
145 -105.00
146 -90.00
147 -75.00
148 -60.00
149 -45.00
150 -30.00
151 -15.00
152 .00
BOUNDARY CONDITIONS
E-15
-------
HEAD IN IMPOUNDMENT H= 100. 00
UNDERLYINB SOIL SUCTION PRESSURE PSIBOT= .000
STEADY STATE ALGORITHM TOOK 23 ITERATIONS, ERROR= 7.3772E-02
TIME(YRS)= .00 TIME STEP= 0 DT=
ITER= 0 TOTAL K ITERATIONS=
PARTICLE DEPTH IN LINER= .000E+00
1.000E+02 ERR"
23
.000E+0B
NODE
POTENTIAL
MOISTURE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
1 . 0000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E-1-02
.4950
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
1 . 0000E-07
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
J . 7230E-09
1 . 7230E-09
1 . 7230E-09
1 . 7230E-09
E-16
-------
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.2682
.2122
.2122
.2122
.2122
.2122
.2682
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7488E-07
7750E-05
7750E-05
7750E-05
7750E-05
7750E-05
7488E-07
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
1.7230E-09
1.7230E-09
E-17
-------
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
-1.
-1
-1
-1
-1
""" I
-1
—
*™
—
1
••»
-1
-1
-1
-1
-1
-1
"""" 1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
-1
— - 1
— 1
— 1
-1
-1
— 1
-1
-1
-1
•H <
-1
-1
-1
-1
-1
-1
-1
-3
-2
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
5000E+02
5000E+02
5000E+02
5000E+02
5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
5000E+02
5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-3.0000E+02
-2.8500E+02
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.2277
.1811
1.
1,
1.
1,
1.
1,
1.
1,
1.
1,
1.
1,
1.
1,
1,
1,
1.
1,
1,
1,
1.
1,
1,
1
1,
1
1,
1
1,
1
1,
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
3
3
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
7230E-09
3.6855E-08
3.0568E-06
E-18
-------
134 -2.7000E+02
135 -2.5500E+02
136 -2.4000E+02
137 -2.2500E+02
138 -2.1000E+02
139 -1.9500E+02
140 -1.8000E+02
141 -1.6500E+02
142 -1.5000E+02
143 -1.3500E+02
144 -1.2000E+02
145 -1.0500E+02
146 -9.0000E+01
147 -7.5000E+01
148 -6.0000E+01
149 -4.5000E+01
150 -3.0000E+01
151 -1.5000E+01
152 .0000E+00
. 1835
. 1861
. 1889
. 1919
.1952
. 1989
.2028
.2072
.2122
.2177
.2242
.2317
.2407
.2518
.2660
.2856
.3157
.3705
.3950
3.5450E-06
4.1461E-06
4.8955E-06
5.8427E-06
7.0588E-06
8.6484E-06
1.0770E-05
1.3670E-05
1.7750E-05
2.3693E-05
3.2719E-05
4.7177E-05
7.1979E-05
1.1863E-04
2.1868E-04
4.8107E-04
1.4615E-03
8.6547E-03
1.7600E-02
TIME(YRS)= .50 TIME STEP= 53 DT=
ITER= 3 TOTAL K ITERATIONS^
PARTICLE DEPTH IN LINER= -1.942E+01
2.820E+05 ERR=
556
6.202E-02
NODE POTENTIAL
1 1.0000E+02
2 9.9762E+01
3 9.9524E+01
4 9.92B7E+01
5 9.9049E+01
6 9.8811E+01
7 9.8573E+01
8 9.8336E+01
9 9.B098E+01
10 9.7860E+01
11 9.7622E+01
12 9.6434E+01
13 9.5245E+01
14 9.4056E+01
15 9.2867E+01
16 9.1679E+01
17 9.0490E+01
18 8.9301E+01
19 B.B112E+01
MOISTURE
K
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
E-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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
8.6924E+01
8.5735E+01
8.3357E+01
8.0980E+01
7.8602E+01
7.6225E-H31
7.3847E+01
7. 1470E+01
6.9092E+01
6.6715E+01
6.4337E+01
6. 1959E+01
5.8393E+01
5.4B27E+01
5.1261E+01
4.7694E+01
4.4128E+01
4.0562E+01
3.6995E+01
3.3429E+01
2.9B63E+01
2.6297E+01
2. 1541E+01
1.67B6E+01
1.2031E+01
7.2763E+00
2.5212E+00
-2.3571E+00
-8. 1784E+00
-1.6487E+01
-3.0123E+01
-5.2751E+01
-8.4043E+01
-1. 1554E+02
-1.390BE+02
-1.5431E+02
-1.6531E+02
-1.6031E+02
-1.5531E+02
-1.5032E+02
-1.4532E+02
-1.4032E+02
-1.3532E+02
-1.3544E+02
-1.3555E+02
-1.3567E+02
-1.357BE+02
-1.3590E+02
-1.3601E+02
-1.3612E+02
-1.3623E+02
-1.3634E+02
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4947
.4855
.4664
.4379
.4020
.3678
.3437
.3298
.3221
.2621
.2087
.2103
.2120
.2138
.2157
.2747
.3318
.3317
.3317
.3316
.3315
.3315
.3314
.3313
.3313
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
9.
7.
4.
2.
1.
4.
2.
1.
1.
1.
1.
1.
1.
1.
2.
2.
2.
2.
2.
2.
2.
2.
2.
2.
2.
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
6468E-0B
5128E-08
4.6619E-08
3107E-08
0029E-0S
6604E-09
7074E-09
9648E-09
6400E-09
4066E-07
4794E-05
6135E-05
764BE-05
9362E-05
1311E-05
2023E-07
0573E-09
2.0543E-09
2.0513E-09
04B3E-09
2.0453E-09
2.0424E-09
2.0394E-09
2.0366E-09
2.0337E-09
E-20
-------
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
-1.3645E+02
-1.3699E+02
-1.3752E+02
-1.3803E+02
-1.3B54E+02
-1.3903E+02
-1.3952E+02
-1.3999E+02
-1.4046E+02
-1.4092E+02
-1.4137E+02
-1.4224E+02
-1.4307E+02
-1.43B5E-M32
-1.4459E+02
-1.4526E+02
-1.45B9E+02
-1.4646E+02
-1.4697E+02
-1.4742E+02
-1.4783E+02
-1.4B34E+02
-1.4875E+02
-1.4907E+02
-1.4932E+02
-1.4951E-M32
-1.4965E+02
-1.4976E+02
-1.4983E+02
-1.49B9E+02
-1.4992E+02
-1.4996E+02
-1.4997E+02
-1.4999E+02
-1.4999E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5000E+02
-1.5001E+02
-1.5002E+02
-1.5004E+02
-1.5007E+02
-1.5012E+02
-1.5021E+02
-1.5036E+02
-1.5061E+02
-1.5101E+02
-1.5164E+02
-1.5261E+02
-1.5408E+02
-1.5624E+02
.3312
.3309
.3306
.3304
.3301
.3298
.3296
.3293
.3291
.3288
.3286
.3281
.3277
.3273
.3269
.3266
.3263
.3260
.3257
.3255
.3253
.3250
.3248
.3247
.3246
.3245
.3244
.3243
.3243
.3243
.3243
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3242
.3241
.3240
.3239
.3237
.3234
.3230
.3223
.3212
2
2
2
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0308E-09
.0170E-09
0036E-09
,9907E-09
9782E-09
,9660E-09
9542E-09
, 9426E-09
9314E-09
, 9205E-09
9099E-09
, 8896E-09
8707E-09
,8530E-09
8367E-09
,B218E-09
8083E-09
,7962E-09
7853E-09
,7757E-09
7673E-09
,7567E-09
7483E-09
.7417E-09
7366E-09
,7328E-09
7300E-09
,7279E-09
7264E-09
,7253E-09
7245E-09
,7239E-09
7235E-09
,7233E-09
7232E-09
.7231E-09
7230E-09
.7230E-09
7229E-09
.7228E-09
,7226E-09
,7223E-09
7217E-09
,7207E-09
7189E-09
.7159E-09
7110E-09
.7031E-09
6907E-09
,6720E-09
6444E-09
1.6050E-09
E-21
-------
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
-1.
-1.
-1.
-1.
-1.
W.O
-2.
-2.
-2.
_o
wO
-2.
-2.
-2,
-2.
-1,
-1,
-1,
-1.
-1
""" 1 »
-1
-8,
-7
-5,
-4
-3
-1
5937E+02
63B3E+02
7007E+02
7B72E+02
-1.9062E+02
2.06B8E+02
-2.2902E+02
-2.590BE+02
-2.9963E+02
.8469E+02
.6975E+02
-2.5480E+02
-2.3984E+02
-2.2487E+02
-2.0990E+02
9493E+02
7994E+02
6496E+02
4997E+02
3498E+02
1999E+02
0499E+02
-8.9996E+01
-7.4998E+01
-5.9999E+01
-4.5000E+01
-3.0000E+01
-1.5000E+01
.0000E+00
.3198
.3178
.3151
.3115
.3069
.3012
.2942
.2860
.2277
.1811
.1835
.1861
.1889
.1920
.1953
. 1989
.2028
.2072
.2122
.2178
.2242
.2317
.2407
.2518
.2660
.2856
.3157
.3705
.3950
1.
1.
1.
1.
1.
9.
8.
6.
3.
3.
3.
4.
4.
5.
7.
8.
1.
1.
1.
2,
3.
4.
7.
1.
2.
4.
1.
8.
1.
5505E-09
4777E-09
3844E-09
2694E-09
1342E-09
8270E-10
8.2217E-10
6.6202E-10
3.6959E-08
3.0659E-06
5541E-06
1553E-06
9047E-06
5.8518E-06
0679E-06
6575E-06
0779E-05
3679E-05
7760E-05
3702E-05
2728E-05
71B6E-05
1988E-05
1864E-04
1868E-04
4.8108E-04
1.4615E-03
8.6547E-03
1.7600E-02
VOLUME BALANCE CALCULATIONS
TOP FLUX 3.3775E-07
BOTTOM FLUX 1.2555E-0B
STORAGE RATE 3.1B35E-07
ERROR
6.B494E-09
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
9.5256E-02
3.5408E-03
B.9783E-02
1.9317E-03
CUMULATIVE CHANGES
VOLUME IN (-)
VOLUME OUT
STORAGE
9.3592E+00
3.1B59E-01
B.3925E+00
E-22
-------
ERROR
RELATIVE ERROR
6.4B16E-01
.077231
TIME(YRS)= 1.00 TIME STEP= 80 DT=
ITER= 1 TOTAL K ITERATIONS=
PARTICLE DEPTH IN LINER= -2.926E+01
3.B00E+04 ERR=
665
5.90BE-02
NODE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
1,
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
9
8
8
8
8
8
8
7
7
7
7
7
6,
6
6,
6
POTENTIAL
0000E+02
9808E+01
9615E+01
9423E+01
9231E+01
9039E+01
S846E+01
B654E+01
8462E+01
8270E+01
B077E+01
7116E+01
6155E+01
5193E+01
4232E+01
3271E+01
2310E+01
134BE+01
0387E+01
8.9426E+01
8.8464E+01
8.6542E+01
8.4619E+01
8.2697E+01
8.0774E+01
8851E+01
6929E+01
5006E+01
3084E+01
1161E+01
9238E+01
6354E+01
3471E+01
6.0587E+01
MOISTURE
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
K
1.0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
1.0000E-07
E-23
-------
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
5.7703E+01
5.4S19E+01
5. 1935E+01
4.9051E+01
4.6167E+01
4.3283E+01
4.0399E+01
3.6554E+01
3.2709E+01
2.8B64E+01
2.501BE+01
2. 1173E+01
1.732SE+01
1.34S3E+01
9.6376E+00
5.7924E+00
1 . 9472E+00
-1.9802E+00
-6.5260E+00
-1.3054E+01
-2.5282E+01
-6.0802E+01
-5.5807E+01
-5.0812E+01
-4.5814E+01
-4.0B16E+01
-3.5817E+01
-3.0818E+01
-3. 1206E+01
-3. 1599E+01
-3. 1998E+01
-3.2402E+01
-3.2812E+01
-3.322BE+01
-3.3650E+01
-3.4077E-M31
-3.4510E+01
-3.4949E+01
-3.7236E+01
-3.96B0E+01
-4.2287E+01
-4.5063E+01
-4.8009E+01
-5. 1127E+01
-5.4414E+01
-5.7864E+01
-6. 1468E+01
-6.5214E+01
-7.3050E+01
-8. 1192E-H31
-8.9414E+01
-9.7469E+01
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4949
.4889
.4743
.4474
.3285
.2708
.2772
.2843
.2926
.3022
.3751
.4359
.4351
.4344
.4337
.4329
.4322
.4314
.4307
.4299
.4291
.4252
.4212
.4170
.4128
.4085
.4042
.3998
.3954
.3910
.3867
.3783
.3704
.3631
.3566
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
9.
8.
5.
2.
1.
2.
3.
4.
6.
8.
5.
2.
2.
2.
2.
2.
2.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
9.
8.
7.
7.
5.
4.
4.
3.
1.0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
7381E-0B
1831E-08
6902E-08
9066E-08
2970E-06
6670E-04
4486E-04
5799E-04
2B57E-04
8.9919E-04
5144E-06
2016E-08
1637E-08
1263E-08
0893E-08
0528E-08
0167E-08
9B10E-08
945BE-08
9110E-08
8766E-0B
7116E-08
557BE-08
4153E-08
2840E-08
1637E-08
0540E-08
5445E-09
8.6454E-09
8372E-09
1135E-09
8955E-09
4.9395E-09
4.1967E-09
3.6240E-09
E-24
-------
87 -1.0512E+02
88 -1.1219E+02
89 -1.1852E+02
90 -1.2407E+02
91 -1.2B82E+02
92 -1.3281E+02
93 -1.3749E+02
94 -1.4092E+02
95 -1.4340E+02
96 -1.451BE-H82
97 -1.4645E+02
98 -1.4737E+02
99 -1.4804E-M32
100 -1.4B52E+02
101 -1.4889E+02
102 -1.4916E+02
103 -1.4942E+02
104 -1.4960E+02
105 -1.4974E+02
106 -1.4984E+02
107 -1.4992E+02
108 -1.5000E+02
109 -1.5008E+02
110 -1.5017E+02
111 -1.5030E+02
112 -1.5046E+02
113 -1.5069E+02
114 -1.5101E+02
115 -1.5146E+02
116 -1.5207E+02
117 -1.52S9E+02
118 -1.5399E+02
119 -1.5545E+02
120 -1.5737E+02
121 -1.5985E+02
122 -1.6303E+02
123 -1.6710E+02
124 -1.7224E+02
125 -1.7S72E+02
126 -1.B6B4E+02
127 -1.9699E+02
128 -2.0966E+02
129 -2.2545E+02
130 -2.4509E+02
131 -2.6949E+02
132 -2.9973E+02
133 -2.8477E+02
134 -2.6982E+02
135 -2.54B5E+02
136 -2.39BBE+02
137 -2.2491E+02
138 -2.0993E+02
.3508
.3459
.3418
.3383
.3355
.3332
.3307
.3288
.3275
.3266
.3260
.3255
.3252
.3249
.3248
.3246
.3245
.3244
.3243
.3243
.3243
.3242
.3242
.3241
.3241
.3240
.3239
.3237
.3235
.3232
.3228
.3223
.3216
.3207
.3196
.3181
.3163
.3142
.3115
.3083
.3046
.3003
.2953
.2896
.2834
.2277
. 1811
. 1835
. 1861
. 1889
. 1920
. 1953
3.
2,
2.
2,
2,
2,
2,
1,
1,
1,
1.
1,
1,
1
1,
1
1,
1
1,
1
1,
1
1,
1
1
1
1,
1
1
1
1
1
1,
1
1
1
1
1
1
1
1
9
8
7
6
3
3
3
4
4
5
7
3.1B46E-09
B4B3E-09
5911E-09
3945E-09
2439E-09
12B3E-09
2.0044E-09
9205E-09
8632E-09
8237E-09
7962E-09
7768E-09
7630E-09
7529E-09
7455E-09
7400E-09
7347E-09
7310E-09
7282E-09
7262E-09
7245E-09
7230E-09
7214E-09
7195E-09
7171E-09
7138E-09
7092E-09
7029E-09
6943E-09
6825E-09
6667E-09
6460E-09
6192E-09
5B50E-09
5424E-09
4903E-09
4277E-09
3541E-09
2695E-09
1746E-09
0707E-09
5993E-10
B.451BE-10
29B6E-10
1768E-10
6931E-08
0634E-06
5516E-06
152BE-06
9022E-06
8494E-06
7.0655E-06
E-25
-------
139 -1.9495E+02
140 -1.7996E+02
141 -1.6497E+02
142 -1.499BE+02
143 -1.3499E+02
144 -1.1999E+02
145 -1.0499E+02
146 -B.9997E+01
147 -7.4999E+01
148 -5.9999E+01
149 -4.5000E+01
150 -3.0000E+01
151 -1.5000E+01
152 .0000E+00
.1989
.2028
.2072
.2122
.2178
.2242
.2317
.2407
.2518
.2660
.2856
.3157
.3705
.3950
8.6551E-06
1.0776E-05
1.3677E-05
1.7757E-05
2.3699E-05
3.2726E-05
4.7184E-05
7.1985E-05
1.1864E-04
2.1868E-04
4.8107E-04
1.4615E-03
S.6547E-03
1.7600E-02
VOLUME BALANCE CALCULATIONS
TOP FLUX 2.9226E-07
BOTTOM FLUX 9.4161E-09
STORAGE RATE 2.8520E-07
ERROR
2.3541E-09
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
1.1106E-02
3.57B1E-04
1.083BE-02
8.945BE-05
CUMULATIVE CHANGES
VOLUME IN (-)
VOLUME OUT
STORAGE
ERROR
RELATIVE ERROR
1.4122E+01
4.8415E-01
1.2812E+01
8.2564E-01
.064443
TIME(YRS)= 2.00 TIME STEP= 133 DT=
ITER= 2 TOTAL K ITERATIONS*
PARTICLE DEPTH IN LINER= -4.744E+01
7.600E+04 ERR=
829
1.153E-02
E-26
-------
NODE
POTENTIAL
MOISTURE
K
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
48
49
1.0000E+02
9.9845E+01
9.9691E+01
9.9536E+01
9.9381E+01
9.9227E+01
9.9072E+01
9.8917E+01
9.8762E+01
9.8608E+01
9.B453E+01
9.7680E+01
9.6906E+01
9.6133E+01
9.5359E+01
9.45B6E+01
9.3812E+01
9.3039E+01
9.2265E+01
9.1492E+01
9.0718E+01
8.9171E+01
8.7624E+01
8.6077E+01
8.4530E+01
B.2983E+01
8.1436E+01
7.9889E+01
7.8342E+01
7.6796E+01
7.5249E+01
7.292BE+01
7.0608E+01
6.B2B7E+01
6.5967E+01
6.3646E+01
6.1326E+01
5.9005E+01
5.6685E+01
5.4364E+01
5.2044E+01
4.8950E+01
4.5856E+01
4.2762E+01
3.9668E+01
3.6574E+01
3.3480E+01
3.0386E+01
2.7293E+01
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
E-27
-------
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
2.'
2. :
1.1
1.'
1.
8.'
5.(
1.1
1.1
2.1
2.!
3.1
3.!
3.1
3.
3.!
3.
3.
3.
3.
3.
3.
3.
3.
3.
3.
3.
3.
2.
2.
2.
2.
2.
2.
2.
2.
2.
1.
1.
1.
1.
1.
1.
8.
6.
3.
1.
-7.
-3.
-5.
-8.
"""" 1 •
4199E+01
1105E+01
B011E+01
4917E+01
1823E+01
8.7290E+00
6350E+00
0635E+01
5635E+01
0635E+01
5635E+01
0635E+01
3.5635E+01
3.5480E+01
5325E+01
5171E+01
5016E+01
4861E+01
4706E+01
4552E+01
4397E+01
4242E+01
4088E+01
3314E+01
2541E+01
1767E+01
0994E+01
0220E+01
9447E+01
B673E+01
2.7900E+01
2.7126E+01
6353E+01
4806E+01
3259E+01
1712E+01
0165E+01
B618E+01
7071E+01
5524E+01
3977E+01
2430E+01
0BB4E+01
8.5631E+00
2427E+00
9223E+00
6019E+00
1853E-01
1545E+00
-5.8554E+00
-B.9565E+00
2653E+01
.4950
.4950
.4950
.4950
.4950
.4950
.4450
.3950
.3950
.3950
.3950
.3950
.4450
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4941
.4902
.4838
.4753
1.1
l.(
1.1
1.1
1.1
1.1
4.
1.'
1.
1.'
1.
1.'
4.
1.1
l.i
1.1
1.
1.1
1.
l.i
1.
l.i
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
1.
9.
8.
7.
5.
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
1952E-05
7600E-02
7600E-02
7600E-02
7600E-02
7600E-02
1952E-05
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
, 0000E-07
0000E-07
,0000E-07
,0000E-07
.0000E-07
,0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.0000E-07
.4222E-0B
8.4513E-08
2003E-08
5.8249E-08
E-28
-------
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
-1.7205E+01
-2.5232E+01
-3.6299E+01
-5. 1053E+01
-6.9112E+01
-8.B577E+01
-1.0676E+02
-1.2167E+02
-1.3275E+02
-1.4053E+02
- .4591E+02
- .4973E-M32
- .5265E+02
- .5510E+02
- .5741E+02
-1.597BE+02
-1.6237E+02
-1.6531E+02
-1.6B70E+02
-1.7266E+02
-1.7729E+02
-1.8272E+02
-1.890BE+02
-1.9654E+02
-2.0528E+02
-2. 1552E+02
-2.2752E+02
-2.4157E+02
-2.5B02E+02
-2.7728E+02
-2.9980E+02
-2.8483E+02
-2.6986E+02
-2.5489E+02
-2.3991E+02
-2.2493E+02
-2.0995E+02
-1.9496E+02
-1.7997E+02
-1.649BE+02
-1.4998E+02
-1.3499E+02
-1.1999E+02
-1.0500E+02
-8.999BE+01
-7.4999E+01
-6.0000E+01
-4.5000E+01
-3.0000E+01
-1.5000E+01
. 0000E+00
.4648
.4475
.4268
.4043
.3824
.3638
.3497
.3398
.3333
.3290
.3263
.3244
.3229
.3218
.3207
.3196
.3184
.3171
.3156
.3140
.3121
.3099
.3075
.3048
.3017
.2984
.2947
.2906
.2862
.2816
.2277
.1811
. 1835
. 1B61
. 1889
.1920
. 1953
.1989
.2028
.2072
.2122
.2178
.2242
.2317
.2407
.2518
.2660
.2856
.3157
.3705
.3950
4
2
1
1,
6
4
3
2,
2
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
9
9
8
7
6
5
3
3
3
4
4
5
7
8
1
1
1
2
3
4
7
1
2
4
1
a
i
4748E-0B
9139E-0B
7765E-0B
0564E-08
4637E-09
2642E-09
1014E-09
4767E-09
1300E-09
9297E-09
B078E-09
7283E-09
6713E-09
6255E-09
5B43E-09
5435E-09
5009E-09
4547E-09
4040E-09
3483E-09
2B74E-09
2213E-09
1503E-09
0750E-09
9610E-10
1462E-10
B.3170E-10
4864E-10
6678E-10
B748E-10
6911E-08
0617E-06
5499E-06
1511E-06
9004E-06
8476E-06
0637E-06
8.6533E-06
1.0775E-05
1.3675E-05
7755E-05
3697E-05
2724E-05
71S2E-05
19B4E-05
1864E-04
1868E-04
8107E-04
4615E-03
B.6547E-03
1.7600E-02
E-29
-------
VOLUME BALANCE CALCULATIONS
TOP FLUX 2.5469E-07
BOTTOM FLUX 7.0621E-09
STORAGE RATE 2.4756E-07
ERROR
7.6461E-11
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
1.9357E-02
5.3672E-04
1.S814E-02
5.B108E-06
CUMULATIVE CHANGES
VOLUME IN (-)
VOLUME OUT
STORAGE
ERROR
RELATIVE ERROR
2.3130E+01
7.2775E-01
2.1579E+01
8.2387E-01
.038180
TIME(YRS)= 3.00 TIME STEP= 186 DT=
ITER= 1 TOTAL K ITERATIONS=
PARTICLE DEPTH IN LINER= -6.215E+01
7.600E+04 ERR*
1001
8.284E-02
NODE POTENTIAL
1 1.0000E+02
2 9.9892E+01
3 9.97B3E+01
4 9.9675E+01
5 9.9567E+01
6 9.9458E+01
7 9.9350E-H31
8 9.9241E+01
9 9.9133E+01
10 9.9025E+01
11 9.8916E+01
12 9.8374E+01
MOISTURE
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1,
1
1
1
1
1
1
1
1
1
1
1
K
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
1.0000E-07
0000E-07
0000E-07
0000E-07
1.0000E-07
E-30
-------
13
14
15
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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
9
9
9
9
9
9
9
9
9
9
9
9
8
8
8
8
8
8
8
8
7
7
7
7
7
7
6
6
6
6
6
5
5,
5
5,
5
4,
4
4,
4
4,
3
3,
3
3,
4
4,
5
5,
6
6.
6,
9.7833E+01
7291E+01
6749E+01
6207E+01
5665E+01
5123E+01
4582E+01
9.4040E+01
3498E+01
2414E+01
1331E+01
0247E+01
9163E+01
B.B0B0E+01
8.6996E+01
8.5912E-I-01
8.4829E+01
B.3745E+01
8.2661E+01
8.1036E+01
9410E+01
7785E+01
6159E+01
4534E+01
290BE+01
12B3E+01
9657E+01
8032E+01
6406E+01
4239E+01
2072E+01
9904E+01
7737E+01
5569E+01
3402E+01
1235E+01
9067E+01
6900E+01
4733E-+-01
2565E+01
0398E+01
8231E+01
6063E+01
3896E+01
8896E+01
3896E+01
8896E+01
3B96E+01
8896E+01
3896E+01
6.3787E+01
3679E+01
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4450
.3950
.3950
.3950
.3950
.3950
.4450
.4950
.4950
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
4. 1952E-05
1 . 7600E-02
1 . 7600E-02
1 . 7600E-02
1 . 7600E-02
1 . 7600E-02
4. 1952E-05
1 . 0000E-07
1 . 0000E-07
E-31
-------
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
5.
4.
4.
4.
4.
4.
4.
4.
4.
3.
3.
3.
3.
3.
3.
2.
2.
2.
2.
1.
1.
1.
1.
1.
8.
6.
4.
2.
-3.
3571E+01
3462E+01
3354E+01
3246E+01
3137E+01
3029E+01
2920E+01
2812E+01
2270E+01
1729E+01
6.1187E+01
6.0645E+01
6.0103E+01
9561E+01
5.9020E+01
847BE+01
7936E+01
7394E+01
6311E+01
5.5227E+01
4143E+01
5.3060E+01
1976E+01
0893E+01
9B09E+01
B725E+01
7642E+01
655BE+01
4933E+01
3307E+01
1682E+01
0057E+01
8431E+01
6B06E+01
51B1E+01
3555E+01
1930E+01
0304E+01
8137E+01
5970E+01
3803E+S1
1636E+01
9469E+01
7301E+01
5134E+01
2967E+01
0800E+01
8.6326E+!Z)0
4654E+00
29B3E+00
1311E+00
-3.6094E-02
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
E-32
-------
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
-2
-4
-7
-1
-1
-2
-3
-4
-7,
-1
-1,
"~ 1
-2,
-2
-2,
-2,
-2,
-2,
-2.
wO
-2.
-2,
-1.
™* 1 i
-1.
-1,
-1.
""" 1 i
w J ^
-8,
-7.
-5,
-4.
-3,
™" 1 •
i
-2.2758E+00
-4.750BE+00
6574E+00
1286E+01
6099E+01
-2.2B60E+01
-3.2825E+01
-4.7880E+01
-7.0218E+01
0086E+02
3767E+02
7574E+02
1080E+02
-2.4176E+02
-2.7027E+02
-2.9974E+02
-2.S479E+02
-2.69B3E+02
-2.5486E+02
3989E+02
-2.2491E+02
-2.0993E+02
9495E+02
7996E+02
6497E+02
499BE+02
3499E+02
1999E+02
0500E+02
-8.9997E+01
-7.4999E+01
-5.9999E+01
-4.5000E+01
0000E+01
5000E+01
0000E+00
.4948
.4921
.4866
.4785
.4673
.4524
.4329
.4087
.3813
.3540
.3306
.3127
.2999
.2906
.2832
.2277
. 1811
. 1835
. 1861
. 1889
. 1920
. 1953
. 1989
.2028
.2072
.2122
.2178
.2242
.2317
.2407
.2518
.2660
.2856
.3157
.3705
.3950
9
B
7
6
4
3
2
1
6
3
1
1
9,
7
6,
3
3,
3,
4,
4,
5,
7,
8.
1,
1.
1,
2.
3,
4.
7,
1.
2,
4.
1.
8.
1.
9.6674E-08
8.8764E-08
7241E-0B
3074E-08
7670E-0B
2873E-0B
0517E-08
1686E-0B
2960E-09
4185E-09
9999E-09
3074E-09
50S4E-10
7.4761E-10
6.1454E-10
6927E-0B
0630E-06
5512E-06
4.1524E-06
4.9018E-06
5.8489E-06
0650E-06
8.6546E-06
0776E-05
3676E-05
7757E-05
3699E-05
2725E-05
71B3E-05
19B5E-05
1864E-04
1868E-04
S107E-04
4615E-03
S.6547E-03
1.7600E-02
VOLUME BALANCE CALCULATIONS
TOP FLUX 2.0837E-07
BOTTOM FLUX 8.6315E-09
STORAGE RATE 1.9916E-07
ERROR
5.8116E-10
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
1.5836E-02
6.5599E-04
1.5136E-02
4.4168E-05
E-33
-------
CUMULATIVE CHANGES
VOLUME IN (-)
VOLUME OUT
STORAGE
ERROR
RELATIVE ERROR
3.0275E+01
9.3159E-01
2.B348E+01
9.9480E-01
.035092
STEADY STATE ACHIEVED DURING TIME STEP 228
YEAR =
3.B1E+00
TIME(YRS)= 3.81 TIME STEP= 228 DT=
ITER= 2 TOTAL K ITERATIONS=
PARTICLE DEPTH IN LINER= -7.483E+01
6.050E+05 ERR=
1127
5.021E-02
NODE POTENTIAL
1 1.0000E+02
2 9.9903E+01
3 9.9806E+01
4 9.9709E+01
5 9.9612E+01
6 9.9515E+01
7 9.9417E+01
8 9.9320E+01
9 9.9223E+01
10 9.9126E+01
11 9.9029E+01
12 9.8544E+01
13 9.B058E+01
14 9.7573E+01
15 9.7087E+01
16 9.6602E+01
17 9.6116E+01
18 9.5631E+01
MOISTURE
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
K
1.0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
1.0000E-07
E-34
-------
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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
9
9
9
9
9
9
9
8
8
8
8
8
B
8
8
8
7
7
7
7
7
7
6
6
6
6
6
6
5
5,
5
5,
5
4
4
4,
4
4,
4
5,
5
6,
6
7,
7
7,
7
7,
7
7,
7
6.
5145E+01
4660E+01
4174E+01
3203E+01
2233E+01
1262E+01
0291E+01
8.9320E+01
8.8349E+01
8.737BE+01
B.6407E+01
8.5436E+01
B.4465E+01
8.3009E+01
8.1552E+01
8.0096E+01
8640E+01
7183E+01
5727E+01
4270E+01
2B14E+01
1358E+01
9901E+01
7959E+01
6017E+01
4076E+01
2134E+01
0192E+01
5.8250E+01
6308E+01
4366E+01
2424E+01
04B2E+01
8541E-M31
6599E+01
4657E+01
2715E+01
0773E+01
5773E+01
0773E+01
5773E-H31
0773E+01
5773E+01
0773E+01
0676E+01
0579E+01
0482E+01
0384E+01
0287E+01
0190E+01
0093E+01
6.9996E+01
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4450
.3950
.3950
.3950
.3950
.3950
.4450
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
4
1
1
1
1
1
4
1
1
1
1
1
1
1
1
1.0000E-07
1.0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
1952E-05
7600E-02
1.7600E-02
1.7600E-02
7600E-02
7600E-02
1952E-05
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
0000E-07
1.0000E-07
E-35
-------
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
6.
5.
5.
5.
5.
5.
5.
5.
5.
4.
4.
4.
4.
4.
4.
4.
3.
3.
3.
3.
3.
2.
2.
2.
2.
2.
1.
1.
1.
1.
1.
9.
7,
5.
3,
1.
9B99E+01
9802E+01
9317E+01
8831E+01
8346E+01
7860E+01
7375E+01
6889E+01
6404E+01
5918E+01
5433E+01
4948E+01
6.3977E+01
3006E+01
2035E+01
1064E+01
0093E+01
9122E+01
8151E+01
5.7180E+01
5.6210E+01
5239E+01
5.3782E+01
2326E+01
0870E+01
9413E+01
7957E+01
6501E+01
5044E+01
3588E+01
2132E+01
0675E+01
8734E+01
6792E+01
4B50E+01
2908E+01
0966E+01
9025E+01
7083E+01
5141E+01
3199E+01
1258E+01
9316E+01
1.7374E+01
5432E+01
3490E+01
1549E+01
6069E+00
6651E+00
5.7233E+00
3.7815E+00
1.8397E+00
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
.4950
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
1 . 0000E-07
E-36
-------
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
-1.0203E-01
-2. 1028E+00
-4.3082E+00
-6.9335E+00
-1.0320E+01
-1.5093E+01
-2.265BE+01
-3.6972E+01
-7.4002E+01
-2.9433E+02
-2.B031E+02
-2.6615E+02
-2.5189E+02
-2.3751E+02
-2. 2304E+02
-2.0848E+02
-1.9385E-H32
-1.7915E+02
-1.6438E+02
-1.4957E+02
-1.3471E+02
-1. 1981E+02
-1.0489E+02
-8.9937E+01
-7.496BE+01
-5.9987E+01
-4.4995E+01
-2.9999E+01
-1.5000E+01
. 0000E+00
.4950
.4948
.4927
.4881
.4807
.4696
.4528
.4257
.3773
.2287
. 1818
.1841
. 1867
.1894
.1924
. 1956
. 1991
.2031
.2074
.2123
.2179
.2243
.2318
. 2407
.2518
.2660
.2856
.3157
.3705
.3950
1.
9.
9.
8.
6.
5.
3,
1.
5,
3.
3,
3,
4,
5.
5,
7,
8
1,
1
1,
2
3,
4
7
1
2
4
1
8
1
0000E-07
7096E-08
0378E-0B
S.0182E-08
66B0E-08
0522E-08
3221E-08
7295E-08
7697E-09
B471E-08
1991E-06
6871E-06
28B1E-06
0374E-06
9B44E-06
2003E-06
8.7897E-06
0911E-05
3B11E-05
7B91E-05
3B33E-05
2859E-05
7316E-05
2117E-05
1877E-04
1881E-04
8120E-04
1.4616E-03
8.6549E-03
1.7600E-02
VOLUME BALANCE CALCULATIONS
TOP FLUX 1.9709E-07
BOTTOM FLUX 1.93B2E-07
STORAGE RATE 3.3105E-09
ERROR
3.6360E-11
VOLUME IN
VOLUME OUT
STORASE VOLUME
ERROR
1.1924E-01
1.1726E-01
2.0028E-03
2.2002E-05
CUMULATIVE CHANGES
VOLUME IN (->
VOLUME OUT
3.5372E+01
3.5412E+00
E-37
-------
STORAGE 3 . 0782E+0 1
ERROR 1 . 04B6E+00
RELATIVE ERROR .034066
#****##*######***###**#****######*#*#*»*#**#*####**####********####»###
BREAKTHROUGH OCCURRED DURING TIME STEP 656 YEAR = 1.20E+01
EXECUTION ENDED NORMALLY - 0 - WARNINGS
E-38
-------
APPENDIX F
SUBROUTINE DESCRIPTION
F-l
-------
This appendix provides a brief description of the subroutines called from
the SOILINER main program. Figure F-l summarizes the order in which
subroutines are called. The main program first calculates the steady state
PSI distribution as a means of comparison to determine the time at which
steady state is achieved during transient simulations. Then, a new PSI
distribution is determined for each time step based on the specified initial
conditions. At the end of each time step the program (1) checks to see if
steady state has been achieved and (2) advances the particle. After steady
state is reached, the program continues to advance the particle until
breakthrough occurs, at which time the program terminates. The following
subroutine descriptions are presented in the order shown in Figure F-l.
FLOPEN - This subroutine is called for the purpose of opening files and
assigning the unit numbers to be associated with the desired input and output
files.
INPUT - All input is read through this subroutine including control, grid
design, soil properties, and initial conditions data. Subroutine NEWBC is
called from INPUT to set the pressure boundary conditions. All input data is
"echoed" to the general output file if LPRINT set equal to T.
CLEAR - The first two calls to this subroutine fill arrays OLDC and VOLD
(the storage array and velocity vector at tn~ , respectively) with zeros.
The third call also zeros the C (storage vector at tn) array, which enables
SOILINER to utilize the same algorithm for both the transient and steady state
solution when solving the set of simultaneous liner equations for the unknown
DPSI. The difference between the two solution strategies is that for steady
state there does not exist either a previous time step or a storage term.
Thus OLDC, VOLDC, and C are eliminated by setting each array to zero. Note
that the C array is zeroed each iteration since it is calculated in subroutine
SPROP.
SPROP - This subroutine is called a number of times for both tne steady
state and transient solution strategies. It is first utilized in the steady
state algorithm within the k-iteration loop. Here, the unsaturated soil
properties SATK (hydraulic conductivity), RMOIST (moisture content), and C
(moisture capacity) are updated using characteristic curves to reflect the
newly calculated PSI distribution. Note that the first time SPROP is called,
soil properties are updated with respect to the specified initial conditions.
When convergence occurs, SPROP is called to update properties consistent with
the final PSI distribution. Finally, SPROP is called to re-establish initial
soil properties for the transient solution. Subsequent calls to SPROP during
a transient simulation serve the identical function described above for steady
state.
BUILDV - This subroutine computes the velocity vector held in the array
VNEW (see equation 3-10 of Part I). For both the steady state and transient
solutions VNEW is incorporated into the forcing vector, F (see the discussion
on BUILDF below). Also, BUILDV is re-initialized at the beginning of the
transient solution, based on the specified initial conditions of PSI, for
subroutine PREDCT.
BUILDF - The right-hand side (RHS) of equation 3-19 (Part I) constitutes
the forcing vector F, of the matrix equation solved using the Thomas algorithm
for the unknown values DPSI. BUILDF builds the forcing vector which includes
VNEW and all other knowns from the RHS of 3-19.
F-2
-------
FLOP EN
INPUT
CLEAR
CLEAR
Iteration
I
CLEAR
BUILDV
BUILDS
BUILDF
THOMAS
FTJROR
SPROP1-*— '
OUTPUT
SPROP1
VCALC
BUILDV
MASBAL
0
•H
I I
rt
cu
steady-state k-it
•rl
BUILDV
BUILDS
BUILDF
THOMAS
SPROPl-^— '
STEADY
> KZOUT
£ OUTPUT
« MASBAL
MF1JTTM
1 '^ NEWBL.
if ERRMAX met
nr MAYTT
exceeded
*__. pnpv
if ERRMAX met
IT- MAYTT
exceeded
SIMPUN
LOCATE
breakthrough
STEADY STATE
SOLUTION
TRANSIENT
SOLUTION
..
Figure F-l. Sequence of SUBROUTINE calls.
F-3
-------
BUILDS - The LHS of equation 3-19 contains the unknown values DPSI for a
given node i and it's two surrounding nodes i-1 and i+1. The LHS also
contains known coefficients which constitute the stiffness matrix STIFF, a
two-dimensional array created by BUILDS. Values assigned to the elements of
STIFF depend upon the solution strategy being evaluated. For the transient
solution, all diagonals of the tridiagonal matrix being represented by STIFF
contain conductivity and grid geometry data. The center diagonal contains
additional data from the storage array, C, and the time step size. For steady
state, all three diagonals contain only grid and conductivity data since the
storage vector, C, is zero for this solution strategy. Once updated, both
STIFF and F are passed to subroutine THOMAS for each k-iteration shown in
Figure F-l.
THOMAS - The Thomas algorithm solves the tridiagonal matrix equation
containing STIFF, F, and the vector of unknowns, DPSI. THOMAS solves for the
unknown values of DPSI, which are passed back to the main program. Once DPSI
is known, PSI is updated thus completing a k-iteration. Another k-iteration
is initiated unless the error for convergence, ERRMAX, is met (see discussion
on subroutine ERROR below), or the maximum number of specified k-iterations,
MAXIT, is exceeded.
ERROR - As long as MAXIT is not exceeded, SOILINER will continue to
iterate during any given time step until convergence is achieved. The
convergence criterion is tested at the end of each k-iteration by passing the
array DPSI to subroutine ERROR. DPSI contains the nodal values of Ai/^ (see
equation 3-19) which are used to update each corresponding value of i]^. If
the solution is converging at a given time step, the Aij^ values will
continue to decrease. Subroutine ERROR "scans" the DPSI array to find the
largest, absolute value (CHMAX) of DPSI after iteration. If CHMAX is less
than the specified error for convergence (ERRMAX), the iterative procedure is
terminated and the steady state solution is obtained. For transient
simulations, when ERRMAX is achieved at a given time level, a new time step
(DT) is determined and the iterative procedure is initiated again. For either
solution strategy, if the maximum number of iterations specified per time step
(MAXIT) is exceeded, a forced exit from the k-iteration loop occurs.
OUTPUT - Depending on the solution strategy chosen, subroutine OUTPUT may
be called to write the steady state solution or initial conditions and
subsequent output for the specified special output times. In either case,
output consists of: (1) time-level data including the specified time (TIME),
time step (NT), time step size (DT), maximum calculated DPSI (ERR), the number
of iterations for the given time step (ITER), the total number of k-iterations
(NKITER) as of time t, and the particle depth (TZDIST); (2) node data,
including pressure, moisture, and conductivity (PSI, RMOIST, and RK,
respectively); and (3) mass balance calculations.
F-4
-------
VCALC - The transient solution begins when the soil properties are
re-established with respect to the specified initial PSI distribution. After
the call to SPROP, the initial particle velocity is determined by passing soil
properties and grid data to subroutine VCALC. Equations 3-25 and 3-26 of
Part 1 are solved for the element flux and velocity respectively in VCALC.
Data required to solve these equations include the conductivity for the
element in which the particle is located, and ty and 8 values from the two
surrounding nodes. For the first call to VCALC, an initial velocity is
returned to the main program in the array VELO. After the initial velocity
calculation, VCALC will be called at the end of each time step. However, the
array VELO will never contain more than 3 velocity values since it is updated
after each call to SIMPUN.
MASBAL - SOILINER provides a mass balance for moisture based upon the
calculated fluxes entering and leaving the flow domain and the change in
storage. Ideally, flux entering the system should equal the rate of storage
plus flux leaving the system. Any difference in this relationship is the
error. MASBAL calculates the flux and storage rates, their corresponding
volumes (per square unit), and the error associated with each time step (note
that the first call to MASBAL is used to initialize variables for subsequent
mass balance calculations). MASBAL also calculates cumulative changes in flux
and storage.
PREDCT - Step one in the finite-difference approximation of the
transient, unsaturated flow equation (see 2-la, Part I) is to predict the new
pressure distribution at time tn+^ = tn +At, based on the soil properties
at tn. Subroutine PREDCT contains equation 3-13 of Part I with ALPHA set
equal to zero, thus making the prediction fully explicit. PREDCT makes the
initial estimate of the new PSI distribution, which is subsequently refined in
the k-iteration loop containing the THOMAS subroutine. The initial estimate
for the nodal values of ipn+l is based on the existing (or initial) PSI
distribution, velocity vector (VOLD) and storage vector (OLDC), and the time
step At.
COPY - This subroutine is called from PREDCT prior to calculation of
3-13. The existing pressure, velocity, and storage vectors (PSI, VNEW, and C,
respectively) are saved into PSIOLD, VOLD, and OLDC. For the initial time
step the arrays PSI, VNEW and C are assigned values in subroutines INPUT,
BUILDV, and SPROP respectively, all occurring prior to subroutine PREDCT. The
COPY function is required not only for equation 3-13 of PREDCT, but also for
subroutine BUILDF in which all six arrays (PSI, PSIOLD, VOLD, C, and COLD) are
required.
STEADY - As a means of determining the time at which steady state is
achieved during a transient simulation, SOILINER first performs the steady
state algorithm. The initial PSI values are read into the array SSPSI, which
ultimately contains the steady state solution. Both the PSI and SSPSI arrays
are passed from the main program to subroutine STEADY where they are compared
node by node. If the maximum difference in PSI of all nodes is less than
1.0 cm, steady state is assumed to have been achieved.
F-5
-------
PTRACK - This subroutine is called at the end of each time step to: (1)
calculate the particle velocity (subroutine VCALC), (2) determine the distance
traveled by the particle during a given time period (subroutine SIMPUN), (3)
advance the particle, and (4) determine the new particle position (subroutine
LOCATE). Each time step, a new particle velocity is determined. For
Simpson's Rule both SIMPUN and LOCATE are called every other time step.
Subroutine PTRACK also passes the integer flag KODEBT to indicate
breakthrough. Initially, the particle is positioned at the liner surface.
When the particle has migrated beyond a given node (NCLAYN - usually specified
at the liner/underlying site-soil interface), KODEBT is assigned the value of
1 to indicate breakthrough, and the program is terminated.
SIMPUN - The objective of this subroutine is to pass a value ZDIST, which
represents the distance traveled by a particle over the period of integration,
to the main program. SIMPUN provides a numerical approximation of the
integral and is capable of handling unevenly spaced points in time. In order
to employ the numerical integration scheme of SIMPUN, both time and
corresponding velocity data are required (VTIME and VELO, respectively). The
current time step is saved in VTIME before each call to VCALC. Simpson's Rule
requires three velocity-time data points to determine a value of ZDIST. Once
returned to the main program, ZDIST is added to the previously calculated,
total distance traveled by the particle (TZDIST) since initiation of a
particular SOILINER simulation.
LOCATE - Each time ZDIST is calculated, it is necessary to determine the
new particle position with respect to the specified grid design. If the
particle passes from one element to the next over the period of integration,
new element data is required for the flux and velocity calculations of
subroutine VCALC. Subroutine LOCATE "scans" the grid design by comparing
successive node depths to the value TZDIST. The first node encountered having
a depth that exceeds TZDIST becomes the bottom node (NODEV2) defining the
element (NELEMV) of the new particle location; the upper node defining NELEMV
(NODEVl) is then calculated as NODEV2-1. If during the scanning procedure the
specified breakthrough node (NCLAYN) is reached and its depth does not exceed
TZDIST, LOCATE sets the integer flag KODEBT equal to one, indicating
breakthrough.
NEWTIM - After the new particle position is established and all output
requirements are satisfied, mass balance calculations are performed (see
subroutine MASBAL). If the specified simulation period or maximum number of
time steps (ENDTIM and MAXNT, respectively) have not been exceeded, a new time
step is then determined. The existing time step (DT), a change parameter
(CHPARM), maximum allowable time step (DTMAX), and current PSI distribution
are passed to subroutine NEWTIM. The first algorithm employed by NEWTIM is
used to determine the maximum change in PSI at any given node (CHMAX) between
successive time steps. Equation 4-1 of Part II is then solved for the new
F-6
-------
DT. The size of DT is due in part to CHMAX. If CHMAX is large, the new DT
will be small, thus preventing excessively large time steps and potential
model inaccuracies. Conversely, if CHMAX is small, DT will increase thus
maintaining model efficiency. Finally, as steady state is approached, CHMAX
will tend to remain relatively small resulting in an ever increasing DT.
Although changes in PSI between successive time steps are small, too large a
DT can still lead to model inaccuracy. To eliminate the potential for
extremely large time steps, a specified time step size limit (DTMAX) is set.
If DT exceeds DTMAX, DT is set equal to DTMAX within NEWTIM.
F-7
-------
APPENDIX G
SOILINER
FORTRAN SOURCE CODE
G-l
-------
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
SOILINER.FOR MAIN PROGRAM
Finite difference of vertical unsaturated infiltration
Dan Goode
RUBS Johnson
Richard Wozmak
June 1983
GCA/Technoloqy Division, Inc.
213 Burlington Road
Bedford, MA 01730
(617) 275-5444
November 1985
##*****#******#**********»*************************************
COMMON/MASSB/STOR,TOTV1,TOTV2,FLUX10,FLUX20
COMMON/DEVICE/IRD,IPRT
COMMON/FILES/IFGRD,IFSOIL,IFINIT,IFPOUT,IFMOUT,IFLUX,IPDOUT
COMMON/INFO/NUMNP,NUMEL,NPM1,NPM2,NUMEL2,AM1,NX,DTV,
» DEPTH,ENDTIM,DT,DTMAX,ALPHA,PSIBOT,H,NZOUT,SRPARM,
* PSINIT,NCPTS,ERRMAX,MAXIT,CHPARMfMAXNTfSYEAR
COMMON/TIMES/TIME,TIME1,NT,ERR,NOUT,TOUT(10),NOUT1,TOUT1
COMMON/TRACK/NODEV1,NODEV2,NELEMV,KODEBT,NCLAYN f TIMEBT,OLDTD f
* TDIST
COMMON/ERRORS/IERR,IWARN
DIMENSION PSI(200),PSIOLD(200),DZ(200),C(200),
* RKC200),WORK(200),SSPSI(200),
# RMOIST(200),
» STIFF(3,200),DPSI(200),F(200),
* Z(200),SATK(200>,PSICRT(200),POR(200>
DIMENSION VNEW(200),VOLD(200),OLDC(200),DZN(200),STARK(200),
* IZOUTU0) ,ISOIL(200) ,BETA(200)
DIMENSION RKL(400),CL(400),RMSTL(400),OLDML(400)
DIMENSION VTIME(3),VELO(3),ZDIST(3)
REAL SRPSI(200)
CHARACTER*20 FILE1,FILE2,FILE3,FILE4
CHARACTER*72 FILEN
LOGICAL LSTEDY,LPRINT
DATA NMAX/200/
IWARN - 0
NKITER = 0
ITER = 0
Open input and output files
CALL FLOPEN
Read and echo problem parameters and initialize
G-2
-------
CALL INPUT
Build velocity vector
CALL BUILDV(VNEW,SSPSI , STARK, DZN,DZ, BETA, RK>
Build stiffness matrix
CALL BUILDS(STIFF,RK,STARK,C,OLDC,DZN,DZ,BETA>
Build forcing vector — boundary condition terms
CALL BU I LDF ( F , C , OLDC , SSPS I , PS I OLD , VNEW , VOLD >
Solve matrix equation for DPS I using Thomas algorithm
CALL THOMAS ( ST I FF , F , DPS I , WORK , NPM2 )
Update PSI using the relaxation parameter SRPARM
DO 125 I=1,NPM2
SRPSI ( I ) =SRPARM*DPSI ( I )
G-3
-------
SSPSI =SSPSI
-------
C Reinitialize soil properties and output initial conditions
CALL SPROP
C Build velocity vector
CALL BUILDV
-------
CALL BUILDS(STIFF,RK,STARK,C,OLDC,DZN,DZ,BETA)
C Build forcing vector — transient and boundary condition terms
CALL BUILDF(F,C,OLDC,PSI,PSIOLD,VNEW,VOLD)
C Solve matrix equation using Thomas algorithm
C (result is incremental change in PSI, DPSI)
CALL THOMAS(STIFF,F,DPSI,WORK,NPM2)
C Update PSI
DO 25 I=1,NPM2
N=I + 1
PSI(N)=PSI(N)+DPSI(I)
25 CONTINUE
C Check solution for convergence
CALL ERROR < DPSI,NPM2,ERR)
IF (ERR.LT.ERRMAX) GOTO 30
IF (ITER.LE.MAXIT) BOTO 20
C **************************
C End of Main Iteration Loop
C *#*********#******#*#*#***
WRITE(IPRT,2000) NT,ERR
WRITE(#,2000) NT,ERR
IWARN=IWARN+1
30 CONTINUE
IAR6=1
CALL SPROP(PSI,RMOIST,C,
* RK,RMSTL,CL,RKL,STARK,SATK,POR,ISOIL,PSICRT)
C Compare transient PSI with steady state SSPSI
CALL STEADY(PSI,SSPSI,KODESS)
C Track particle and check for breakthrough, steady state
CALL PTRACK(STARK,DZ,RMOIST,PSI,Z,VTIME,VELO,ZDIST,ITER,NKITER1
IF (KODEBT.EQ.1) BO TO 34
IF(KODESS.EQ.1) SO TO 40
C Write to output files
34 CONTINUE
IF (IOUT.EQ.1.AND.NOUT.NE.-1) THEN
NOUT1=NOUT1+1
TOUT1=TOUT(NDUT1)
IF(NOUT1.ST.NOUT.OR.TOUT1.LE.0.D0) NOUT»0
END IF
IF (NOUT.EQ.-1.0R.TIME.6E.ENDTIM.OR.NT.6E.MAXNT) IOUT=1
IF (IOUT.EQ.1) THEN
CALL OUTPUT(PSI,RMOIST,RK,Z,RKL,CL,RMSTL,STARK,POR,
* DZ,ITER,NKITER,TDIST)
G-6
-------
NTOUT=NTOUT+1
END IF
C Calculate mass balance of -flux
CALL MASBAL < PS I , RMSTL , OLDML , STARK , DZ , I OUT )
IF
1900 FORMATC3X, 'T-STEP',4X, 'DT(S) ',5X, 'TIME(Y) ',5X, *ZDIST',6Xf 'TDIST',4
*X, 'NELEMV',2X, 'ITER' ,2X, 'NKITER'/)
1990 FORMAT (//72(1H*)/ IX, 'STEADY STATE ALGORITHM TOOK ',15,' ITERATIONS,
* ERROR= ' ,1PE10.4/72<1H*)//)
2000 FORMAT
-------
c
C SOILINER SUBPROBRAMS
C
C Dan Goode June 1983
C
C BCA/Technology Division, Inc
C 213 Burlington Road
C Bedford, MA 01730
C (617) 275-5444
C
C RUBS Johnson November 1985
C Richard Wozmak
C
£*««***«*********«***««****«*«****«*«*****«*****«******«•************
C
SUBROUTINE FLOPEN
C
C
C Assigns -file unit specifiers (fills DEVICE and FILES commons)
C and opens all i /o files.
C
C******************************************************************
COMMON/DEVICE/ I RD, IPRT
COMMON/FILES/ I F6RD, IFSOIL, IFINIT, IFPOUT, IFMOUT, IFLUX, IPDOUT
CHARACTER*72 FILEN
C Assign file unit specifiers
IRD - 5
IPRT = 7
IFGRD » 10
IFSOIL - 11
IFINIT = 12
IFPOUT - 13
IFMOUT =14
IFLUX • 15
IPDOUT =16
INFN - 17
C Open all output files
OPEN ( I FMOUT , F I LE= ' MO I ST . PRN ' , ST ATUS= ' NEW ' )
OPEN < I FPOUT , F I LE= ' PS I . PRN ' , STATUS= ' NEW ' )
OPEN ( I FLUX , F I LE= ' FLUX . OUT ' , STATUS= ' NEW ' )
OPEN ( I PDOUT , F I LE= ' PDT . OUT ' f STATUS= ' NEW ' )
OPEN ( I PRT , F I LE= ' BEN . OUT ' , ST ATUS= ' NEW ' )
C Open INFN. DAT, read input file names, and open them
OPEN ( INFN,FILE= ' INFN. DAT ' )
READ < INFN, 201) FILEN
G-8
-------
OPEN ( I RD , F I LE=F I LEN , ST ATUS= ' OLD ' )
READ(INFN,201) FILEN
OPEN ( I FBRD , F I LE=F I LEN , ST ATUS= ' OLD ' )
READ (INFN, 201) FILEN
OPEN < I FSO I L , F I LE=F I LEN , ST ATUS= ' OLD ' )
READ (INFN, 201) FILEN
OPEN (IFINIT,FILE=FILEN,STATUS=' OLD')
201 FORMAT TV,
ZOUT,SRPARM,
,SYEAR
TIME=0.D0
TIME1=0.D0
SYEAR=3.1536E+7
NT=0
READ(IRD,628) TITLE
628 FORMAT(A80)
WRITE(IPRT,2001) TITLE
Read temporal control parameters
READ(IRD,*) LPRINT,LSTEDY,ENDTIM,DT,DTMAX,MAXNT,ALPHA,ERRMAX,
* MAXIT,CHPARM
G-9
-------
IF
END IF
Read spatial control parameters
READ
-------
199 CONTINUE
IF (UPLIM.LT.NUMNP) STOP 'LT NUMNP nodes specified in grid file
IF
49 CONTINUE
Echo grid data
IF (LPRINT) THEN
WRITE(IPRT,2091)
WRITE(IPRT,2093)
WRITE(IPRT,2069) 1,Z(1)
WRITE
DO 69 I=2,NPM1
WRITE(IPRT,2069) I,Z(I),DZN(1-1)
WRITE(IPRT,2071) I,DZ(I)
69 CONTINUE
WRITE(IPRT,2069) NUMNP,Z 0'
UPLIM = LOWLIM + NUM - 1
DO 169 K - LOWLIM, UPLIM
ISOIL(K) - ISOILX
SATK(K) = SATKX
POR(K) = PORX
PSICRT(K) = PSICRX
IF (LPRINT) WRITE(IPRT,2070) K,ISOILX,SATKX,PORX,PSICRX
169 CONTINUE
LOWLIM = UPLIM + 1
GOTO 168
299 CONTINUE
IF (UPHM.LT.NUMEL) STOP 'LT NUMEL elements in properties file*
IF (UPLIM.GT.NUMEL) STOP '6T NUMEL elements in properties file'
Read initial conditions
IF (LPRINT) WRITE(IPRT,2089)
LOWLIM = 1
54 CONTINUE
READ(IFINIT,*,END=399> NUM,PSIQ
UPLIM = LOWLIM + NUM - 1
DO 56 J = LOWLIM, UPLIM
PSI(J) = PSIQ
G-ll
-------
IF
1030 FORMAT(15)
1040 FORMAT(3F10.0)
1050 FORMAT(8F10.0)
1061 FORMAT(15,F10.0)
1070 FORMAT(215,3F10.0)
1080 FORMAT(I10,F10.0)
1081 FORMAT(F10.0)
2000 FORMAT(//' TEMPORAL DISCRETIZATION PARAMETERS'//
*' STEADY STATE PARM',12(1H.),'LSTEDY=',L10/
*' IF LSTEDY EQ T, COMPUTE STEADY STATE ONLY'/
*' OTHERWISE, COMPUTE TRANSIENT SOLUTION'/
*' SIMULATION TIME',19(1H.>,'ENDTIM(YRS)*',F10.2/
*' TIME STEP',30(1H.),'DT=',E10.3/
*' MAXIMUM ALLOWABLE TIME STEP',10(1H.>,'DTMAX=',E10.4/
»' MAXIMUM NUMBER OF TIME STEPS',12(1H.),'MAXNT=',I10/
*' TEMPORAL WEIGHTING PARAMETER',10(1H.),'ALPHA*',F10.2/
*' MAXIMUM ERROR FOR CONVERGENCE',7(1H.),'ERRMAX=',F10.4/
*' MAXIMUM ITERATIONS PER TIME STEP',6(1H.>,'MAXIT=*,I10/
*' TIME STEP CHANGE PARAMETER',15(1H.)f'CHPARM=',F10.4)
2001 FORMAT(//IX,80(1H*)//' SOILINER OUTPUT'//
*1X,80(1H*)//1X,A80//1X,80(1H*))
G-12
-------
2004 FORMAT (-) <-) '/>
2040 FORM AT (3F 13. 3)
2050 FORMAT , 'POR= ' ,F10. 4)
2060 FORMAT
C
Computes vector V = C*dPSI/dT
using centered finite difference method
-------
DIMENSION V<1) ,PSI<1) ,STARK<1) ,DZN<1) ,DZ(1) ,RK(1> ,BETA(1)
DO 10 L=1,NPM2
LP1=L+1
V(L)=< STARK /DZN(L)
10 CONTINUE
RETURN
END
C******#***»**************************************************
c
SUBROUT I NE PREDCT < PS I , PS I OLD , VNEW , VOLD , C , OLDC , OLDML , RMSTL )
C
Stores old vectors PSIOLD, VOLD, OLDC, and OLDML,
, , ,
and predicts new PSI values.
***********#*****#******#**#####*****#*******
COMMON/MASSB/STOR , TOTV 1 , TOTV2 , FLUX 1 0 , FLUX20
COMMON / DE V I CE / I RD , I PRT
COMMON / I NFO / NUMNP , NUMEL , NPM 1 , NPM2 , NUMEL2 , AM
# DEPTH , ENDT I M , DT , DTMAX , ALPHA , PS I
* PS I N I T , NCPTS , ERRMAX , MAX I T , CHPAR
, , ,
S I BOT , H , NZOUT , SRPARM ,
ARM , MAXNT , SYEAR
10) ,NOUT1,TOUT1
,,,,
ORS/ I ERR, I WARN
PSI (1) ,PSIOLD<1) ,VNEW(1) ,
OLDML(l) , RMSTL U>
Store old vectors
CALL COPY < PS I , NUMNP , PS I OLD )
CALL COP Y( VNEW, NPM2, VOLD)
CALL COPY < C , NUMNP , OLDC )
CALL COP Y( RMSTL, NUMEL2, OLDML)
Predict new PSI from V
DO 10 I=1,NPM2
N=I + 1
IF (OLDC(N) .BT.0.D0) PSI
-------
c
c***********»***********#**********#**#*********##*#
COMMON/MASSB/STOR,TOTV1,TOTV2,FLUX10,FLUX20
COMMON/DEVICE/IRD,IPRT
COMMON/INFO/NUMNP,NUMEL,NPM1,NPM2,NUMEL2,AM1,N
* DEPTH f ENDTIM,DT,DTMAX,ALPHA,PSIBOT
Loop over interior nodes
DO 10 L=1,NPM2
LP1=L+1
ADZ I NV=ALPHA/DZN < L )
STIFF < 1 f L) =-ADZINV*STARK (L) /DZ (L)
STIFF (2,L) =ADZINV» (STARK (L) /DZ (L) +STARK (LP1 ) /DZ (LP1 > ) +
* ( ALPHA#C ( LP 1 ) + AM 1 *OLDC ( LP 1 ) ) /DT
STIFF <3,L)=-ADZINV*STARK(LP1>/DZ
10 CONTINUE
RETURN
END
C
C
SUBROUT I NE BU I LDF ( F , C , OLDC , PS I , PS I OLD , VNEW , VOLD )
C
C Builds forcing vector F with
C boundary conditions and transients.
C
******»****#*#****************»*#*#**
COMMON/MASSB/STOR , TOTV 1 , TOTV2 , FLUX 1 0 ,
COMMON/DE V I CE / I RD , I PRT
COMMON/ I NFO/NUMNP , NUMEL , NPM 1 , NPM2 , NU
* DEPTH , ENDT I M , DT , DTMAX , AL
DO 10 L=1,NPM2
LP1=L+1
F(L)=ALPHA*VNEW(L) + AM1*VOLD(L) -
* +
* AM1*OLDC(LP1) )/DT
10 CONTINUE
RETURN
END
G-15
-------
c
SUBROUTINE NEWT IM(RNEW, OLD, KODE)
C
C
C Computes new time step DT based on maximum change of PS I
C during previous step, CHMAX.
C
C*******************************#************************#*
COMMON / DE V I CE / I RD , I PRT
COMMON/TIMES/TIME, TIME1, NT, ERR, NOUT, TOUT (10) ,NOUT1,TOUT1
COMMON/ I NFO/NUMNP , NUMEL , NPM 1 , NPM2 , NUMEL2 , AM 1 , NX , DTV ,
* DEPTH , ENDT I M , DT , DTMAX , ALPHA , PS I EOT , H , NZOUT , SRPARM ,
* PS I N I T , NCPTS , ERRMAX , MAX I T , CHPARM , MAXNT , SYEAR
COMMON /ERRORS/ I ERR , I WARN
DIMENSION RNEW(l) ,QLD(1)
C Assign absolute value of largest change in PSI to CHMAX
CHMAX=0.D0
DO 10 I =2, NUMEL
CDP CHANGE NEXT CARD FOR DOUBLE/SINGLE PRECISION
CHANGE-ABS ( RNEW ( I ) -OLD < I ) )
IF < CHANGE. GT. CHMAX) CHMAX =CHANGE
10 CONTINUE
C Compute new time step DT as product of old time step
C times ratio of allowed change to actual change at last step
DT=DT#CHPARM/CHMAX
IF(DT.GT.DTMAX) DT=DTMAX
IF(NOUT.LE.0) GOTO 20
TIM=TIME1+DT
IF (TIM.LT.TOUT1) GOTO 20
DT=TOUT1-TIME1
KODE=1
20 CONTINUE
TIM=TIME1+DT
IF(TIM.GT.ENDTIM) DT=ENDTIM-TIME1
RETURN
END
C
SUBROUT I NE OUTPUT (PS I , RMO I ST , RK , Z , RKL , CL , RMSTL , STARK , POR , DZ ,
* ITER,NKITER,TDIST)
C
C**********************************************»***************
C
C Write time stepping, pressure, moisture, and conductivity data
C to general output file. Also write pressure, moisture, and
C flux data to corresponding output files.
C
C*#***********#***********»**********#**#*****#****##******#**»
G-16
-------
COMMON/DEVICE/IRD,IPRT
COMMON/TIMES/TIME,TIMEl,NT,ERR,NOUT,TOUT(10),NOUTl,TOUT1
COMMON/INFO/NUMNP,NUMEL , NPM1,NPM2,NUMEL2,AMI,NX,DTV,
* DEPTH,ENDTIM,DT,DTMAX,ALPHA,PSIEOT,H,NZOUT,SF
* PSINIT,NCPTS,ERRMAX,MAXIT,CHPARM,MAXNT,SYEAR
COMMON/ERRORS/IERR,IWARN
COMMON/FILES/IFORD,IFSOIL,IFINIT,IFPOUT,IFMOUT,IFLUX,IPDC
Write assorted info, to general output file (GEN. OUT)
YTIME=TIME/SYEAR
WRITE (IPRT, 1900) YTIME,NT,DT,ERR,ITER,NKITER,TDTST
WRITE (IPRT, 2010) (I,PSI(I) ,RMOIST(I> ,RK(I) ,1=1, NUMNP)
C Write pressure data to pressure output file (PSI.PRN)
WRITE (IFPOUT, ' (/) ')
DO 4 1=1, NUMNP
CDP CHANGE NEXT CARD FOR DOUBLE/SINGLE PRECISION
IF (PSI(I) .LT.0.D0) THEN
PF=ALOG10(-PSI(I) )
ELSE
PF=0.D0
END IF
WRITE (IFPOUT, 2040) I,Z(I) ,PSI (I) ,PF
4 CONTINUE
C Write moisture data to moisture output file (MOIST. PRN)
WRITE (IFMOUT, '(/)')
WRITE (IFMOUT, 2030) (L,Z(L) ,RMOIST(L) ,L=1, NUMNP)
C Calculate flux and particle velocity and
C write to flux output file (FLUX. OUT).
DO 40 L=l, NUMEL
LP1=L+1
FLUX=-STARK (L) * ( 1 . D0+ (PSI (LP1 ) -PSI (L> ) /DZ (L) )
VEL=FLUX/ ( (RMOIST (L) +RMOIST (LP1 ) ) /2. D0)
ZZ=(Z(L)+Z(LP1))/2.D0
WRITE (I FLUX, 2040) L, ZZ ,FLUX , VEL
40 CONTINUE
RETURN
1900 FORMAT (//1X, 71 (1H*)//' TIME (YRS) = ' ,F10. 2, ' TIME STEP= ' , 15, ' DT= '
*,1PE12.3,' ERR=',1PE12.3//' ITER= ', 15, 10X, 'TOTAL K ITERATIONS*' , I
*10//' PARTICLE DEPTH IN LINER= ' , 1PE12.3//1X,71 ( 1H*) //
*' NODE POTENTIAL MOISTURE K'/)
2010 FORMAT(I6,4X,1PE13.4,4X,0PF12.4,4X,1PE13.4)
2030 FORMAT(I10,F10.3,F10.7)
2040 FORMAT(I10,0PF10.2,2(1PE10.2,1X,1PE10.3))
END
G-17
-------
c
SUBROUTINE ERROR < A, N, ERRM)
C
C*******************************************************************
C
C Compute scalar error ERRM from maximum absolute value of
C first N components of vector A.
C
DIMENSION Ad)
ERRM=0.D0
DO 10 1=1, N
AERR=ABS
DO 10 1=2, N
10 CONTINUE
C Substitution
C < 1 ) =B < 1 )
DO 30 1=2, N
C ( I ) =B ( I ) -A < 1 , 1 ) *C ( IM1 ) /W ( IM1 )
30 CONTINUE
C(N)=C(N)/W
-------
c
SUBROUT I NE MASBAL ( PS I , RMSTL , OLDML , STARK , DZ , I OUT >
C
C***********************************************^
C
C Computes and prints mass balance of flux calculations.
C
COMMON/MASSB/STOR , TOTV 1 , TOTV2 , FLUX 1 0 , FLUX20
COMMON/TIMES/TIME,TIME1,NT,ERR,NOUT,TOUT(10) ,NOUT1,TOUT1
COMMON/DEVICE/ IRD t IPRT
COMMON/ I NFO/NUMNP , NUMEL , NPMl , NPM2 , NUMEL2 , AMI , NX , DTV ,
* DEPTH,ENDTIM,DT,DTMAX,ALPHA,PSIBOT,H,NZOUT,SRPARM,
* PS I N I T , NCPTS , ERRMAX , MAX I T , CHPARM , MAXNT , SYEAR
COMMON/ERRORS/ I ERR , I WARN
DIMENSION PSI (1) f RMSTL (1) , STARK (1) ,OLDML(1) ,DZ(1)
C Initialization
IF(IOUT.EQ.-l) THEN
TOTV1=0.D0
TOTV2=0.D0
STOR=0.D0
FLUX 10=STARK(1)*(1.D0+ (PSI (2) -PSI (1))/DZ(1) )
FLUX20=STARK (NUMEL) * ( 1 . D0+ (PSI (NUMNP) -PSI (NUMEL) ) /DZ (NUMEL) )
RETURN
END IF
C Top -flux
FLUX 1=STARK (!)*(!. D0+ (PSI (2) -PSI ( 1 ) ) /DZ ( 1 ) )
VOL 1 = ( FLUX 1 *ALPHA+FLUX 1 0*AM 1 ) *DT
FLUX10=FLUX1
TOTV1=TOTVH-VOL1
C Bottom flux
FLUX2=STARK (NUMEL) * ( 1 . D0+ (PSI (NUMNP) -PSI (NUMEL) > /DZ (NUMEL) )
VOL2= (FLUX2*ALPHA+AM1*FLUX20) *DT
FLUX20=FLUX2
TOTV2=TOTV2+VOL2
C Change in storage (assumes linear variation)
DSTOR=0. D0
DO 10 L=l, NUMEL
CDP CHANBE NEXT CARD FOR DOUBLE/SINGLE PRECISION
DELZ=ABS(DZ(L) )
L2 = L*2
LI « L2-1
DSTOR=DSTOR+ (RMSTL (LI ) -OLDML (LI ) +RMSTL (L2) -OLDML (L2) ) *DELZ/
* 2.D0
10 CONTINUE
STOR*STOR+DSTOR
IF(IOUT.NE.l) RETURN
RATEST=DSTOR/DT
G-19
-------
C Compute flux and volume errors EFLUX and EVOL
CDP CHANGE NEXT THREE CARDS FOR DOUBLE/SINGLE PRECISION
EFLUX=ABS(-FLUX1+FLUX2+RATEST)
EVOL=ABS < DSTOR-VOL1+VOL2)
ETOT=ABS < STOR-TOTV1+TOTV2)
EREL=ETOT/STOR
WRITE FLUX1,FLUX2,RATEST,EFLUX,VOL1,VOL2,
* DSTOR,EVOL,TOTV11TOTV2,STOR,ETOT,EREL
RETURN
2020 FORMAT (//' VOLUME BALANCE CALCULATIONS'//
*' TOP FLUX ',1PE12.4/
*' BOTTOM FLUX' ,1PE12.4/
*' STORAGE RATE',1PE12.4/
*20X,12(1H-)/
*' ERROR ',1PE12.4///
*' VOLUME IN ',1PE12.4/
*' VOLUME OUT ',1PE12.4/
* ' STORAGE VOLUME ' , 1 PE 1 2 . 4 /
*20X,12(1H-)/
*' ERROR ',1PE12.4///
*' CUMULATIVE CHANGES'//
#' VOLUME IN (-) ',1PE12.4/
*' VOLUME OUT ',1PE12.4/
*' STORAGE ',1PE12.4/
*20X,12(1H-)/
*' ERROR ',1PE12.4//
*' RELATIVE ERROR ',0PF12.6>
END
C*****************************************************************
C
SUBROUTINE NEWBC(PSI)
C
C
C Reads pressure boundary conditions PSI(l) and PSI (NUMNP)
C from control file.
C
C**************************************************^
COMMON/ 1 NFO/NUMNP , NUMEL , NPM1 , NPM2 , NUMEL2 , AMI , NX , DTV f
# DEPTH , ENDT I M , DT , DTMAX , ALPHA , PS I BOT , H , NZOUT , SRPARM ,
* PSINIT,NCPTS,ERRMAX,MAXIT,CHPARM,MAXNTfSYEAR
COMMON/DEVICE/ I RD, I PRT
DIMENSION PSI (1)
READdRD,*) H7PSIBOT
WRITE (IPRT, 2020) HfPSIBOT
PSK1)=H
PSI
-------
2020 FORMAT
-------
C Evaluate second point (velocity as a function of time) to be
C used in Simpson's method
IF NT,TYRS,TDIST
WRITE<*,1001) NT,DT,TYRS,ZDIST<3),TDIST,NELEMV,ITER,NKITER
IF(KODEBT.EQ.1> RETURN
C Reassign third point of previous calculation to first point
C of subsequent calculation
NODEV1=NELEMV
NODEV2=NODEV1+1
VELO(1)=VELO<3)
VTIME<1)=VTIME<3)
RETURN
900 FORMATC/3X,'T-STEP',4X,'DT(S)',5X,'TIME(Y)',5X,'ZDIST',6X,'TDIST',
*4X,'NELEMV,2X,'ITER',2X,'NKITER'/)
1000 FORMAT(3<5X,1PE10.3>>
1001 FORMAT<4X,15,4<1X,1PE10.3),2X,15,2X,15,3X,15)
1002 FORMAT(I10,2(5X,1PE10.3))
END
C
SUBROUTINE VCALC(VELO,STARK,DZ,RMOIST,PSIf M)
C
C
C Calculates particle velocity, VELO(M), where M corresponds
C to the point number in Simpson's integration scheme.
C
COMMON/TRACK/NODEVl,NODEV2,NELEMV,KODEBT,NCLAYN,TIMEET,OLDTD,
* TDIST
DIMENSION STARK<1>,DZ(1),RMOIST<2),PSI(2),VELO(1)
G-22
-------
FLUX=-STARK(NELEMV)*(1.D0+(PSI(NODEV2)-PSI(NODEV1))/DZ(NELEMV))
VELO(M)=FLUX/((RMOIST(NODEV1)+RMOIST(NODEV2)>/2.0)
RETURN
END
C
C
SUBROUTINE SIMPUN(XX,FX,AX)
C
C
C
C Numerically evaluates the velocity integral (Simpson's rule),
C
C SUBPROGRAM AUTHOR: J. BARISH, COMPUTING TECHNOLOGY CENTER,
C UNION CARBIDE CORP., NUCLEAR DIVISION, OAK RIDGE TENN.
C
C ****************«***********************«^******^
COMMON/INFO/NUMNP,NUMEL,NPM1,NPM2,NUMEL2,AM1,NX,DTV,
* DEPTH,ENDTIM,DT,DTMAX,ALPHA,PSIEOT,H,NZOUT,SRPARM,
* PSINIT,NCPTS,ERRMAX,MAXIT,CHPARM,MAXNT,SYEAR
DIMENSION XX(2), FX(2), AX(2)
AX(l)-0.0
DO 10 IX=2,NX,2
D1=XX(IX)-XX(IX-1)
AX(IX)=AX(IX-1)+D1/2.0*(FX(IX)+FX(IX-1))
IF(NX.EQ.IX) RETURN
D2=XX(IX+1)-XX(IX-1)
D3=D2/D1
A2=D3/6.0*D2#*2/(XX(IX+l)-XX(IX))
A3=*D2/2. 0-A2/D3
10 AX(IX+l)=AX(IX-1)+(D2-A2-A3)»FX(IX-1)+A2*FX(IX)+A3*FX(IX+l)
RETURN
END
C
C
SUBROUTINE LOCATE(Z,VELO,VTIME)
C
C
C
C Updates particle position in grid NELEMV, and
C sets breakthrough flag KODEBT = 1 when particle passes
C through liner base C TDIST > Z(NCLAYN) 3.
C
COMMON/DEVICE/IRD,IPRT
COMMON/FILES/IFGRD,IFSOIL,IFINIT,IFPOUT,IFMOUT,IFLUX,IPDOUT
COMMON/TIMES/TIME,TIME1,NT,ERR,NOUT,TOUT(10),NOUT1,TOUT1
COMMON/TRACK/NODEV1,NODEV2,NELEMV,KODEBT,NCLAYN,TIMEET,OLDTD,
* TDIST
COMMON/INFO/NUMNP,NUMEL,NPM1,NPM2,NUMEL2,AM1,NX,DTV,
* DEPTH,ENDTIM,DT,DTMAX,ALPHA,PSIEOT,H,NZOUT,SRPARM,
G-23
-------
* PSINIT,NCPTS,ERRMAX,MAXIT,CHPARM,MAXNT,SYEAR
DIMENSION Z<1) ,VELO<3) ,vTIMEU>
DATA AVEL/0.0/,CDIST/0.0/
IF*1.175D6/(1.175D6+PSINEG**4)
CL(IP) = 6.3796D6*(POR(L)-0.075D0)*PSINEG**2.96D0/
* ((1.61D6+PSINE6**3.96D0)**2)
SOTO 20
C Soil type 1
101 CONTINUE
RMSTL(IP) - POR(L)/((PSINEB/12.1)**0.2469)
G-26
-------
RKL(IP) = SATK(L)*(RMSTL(IP)/POR(L))**!!.10
102
103
104
105
106
107
108
109
CL(IP) = ,
GOTO 20
Soil type 2
CONTINUE
RMSTL(IP)
RKL(IP) =
CL(IP) *
GOTO 20
Soil type 3
CONTINUE
RMSTL(IP)
RKL(IP) -
CL(IP) =
GOTO 20
Soil type 4
CONTINUE
RMSTL(IP)
RKL(IP) =
CL(IP) -
GOTO 20
Soil type 5
CONTINUE
RMSTL(IP)
RKL(IP) =
CL(IP) »
GOTO 20
Soil type 6
CONTINUE
RMSTL(IP)
RKL(IP) =
CL(IP) =
GOTO 20
Soil type 7
CONTINUE
RMSTL(IP)
RKL(IP) =
CL(IP) =
GOTO 20
Soi1 type 8
CONTINUE
RMSTL(IP)
RKL(IP) «
CL(IP) =
GOTO 20
Soil type 9
CONTINUE
RMSTL(IP)
RKL(IP) =
CL(IP) =
GOTO 20
Soil type 10
4570*POR(L)/PSINEG**!.2470
= POR(L)/((PSINEG/ 9.0)**0.2283)
SATK(L)*(RMSTL(IP)/POR(L))**11.76
3770*POR(L)/PSINEG**1.2283
= POR(L)/((PSINE6/21.8)**0.2041)
SATK(L)*(RMSTL(IP)/POR(L))**12.80
1088*POR(L)/PSINEG**!.2041
= POR(L)/((PSINEG/78.6)**0.18B7)
SATK(L)*(RMSTL(IP)/POR(L))**13.60
4299*POR(L)/PSINEG**!.1887
- POR(L)/((PSINEG/47.B)**0.1855)
SATK(L)*(RMSTL(IP)/POR(L))**13.78
3802*POR(L)/PSINEG**!.1855
= POR(L)/((PSINEG/29.9)**0.1405)
SATK(L> *(RMSTL(IP)/POR(L))**17.24
2264*POR(L)/PSINEG**!.1405
- POR(L)/((PSINEG/35.6)**0.1290)
SATK(L)*(RMSTL(IP)/POR(L))**18.50
2046*POR(L)/PSINEG**1.1290
- POR(L)/((PSINEB/63.0)**0.1174)
SATK(L)*(RMSTL(IP)/POR(L))**20.04
1909*POR(L)/PSINEG**!.1174
= POR(L)/((PSINEG/15.3)**0.0962)
SATK(L)*(RMSTL(IP)/POR(L))**23.80
1250*POR(L)/PSINEG**1.0962
G-27
-------
110 CONTINUE
RMSTL(IP) » POR(L)/((PSINEB/49.0)#*0.0962)
RKL(IP) - SATK(D* (RMSTL(IP) /POR(L) )**23.80
CL(IP) = .1398*POR(L)/PSINEG**1.0962
GOTO 20
Soil type 11
111 CONTINUE
RMSTL(IP) = POR(L)/((PSINEG/40.5)#*0.0877>
RKL(IP) = SATK(L)*(RMSTL(IP)/POR(L))**25.80
CL(IP) = .1214*POR(L)/PSINEB**1.0877
GOTO 20
Soil type 12
112 CONTINUE
RMSTL*124.6D0/(124.6D0+PSINES**!.77D0)
CL(IP) = 2956.D0*(POR(L)-.124D0)*PSILN**3/PSINEB/
# (739.D0+PSILN**4)**2
ELSE
Point of saturation
IF (PSI(N).GE.0.0) THEN
RMSTL(IP) = POR(L)
CL(IP) = 0.0D0
RKL(IP) = SATK(L)
GOTO 20
END IF
Relationship for PSI > PSI at inflection point
PSINEG = -PSKN)
Branch on soil type
BOTO (201,202,203,204,205,206,207,208,209,210,
* 211,212) JSOIL
IF (JSOIL.NE.0) STOP 'JSOIL out of bounds in SPROP'
Soil type 0
RMSTL(IP) » POR(L)
RKL(IP) = SATK(L)
CL(IP) = 0.0D0
GOTO 20
Soil type 1
201 CONTINUE
RMSTL (IP) = POR (L>*. 89850+ (POR (D* (9517112.6-
* (9395661.6+(6867.*PSINE6)))**.5)/3433.S
RKL(IP) = SATK (D* (RMSTL (IP) /POR (L) )**!!. 10
CL(IP) = POR(L>/(9517112.6-
* (9395661.6+(6867.0*PSINE6)))**0.5
GOTO 20
Soil type 2
202 CONTINUE
G-28
-------
RMSTL(IP) = POR (L>». 89540+ (POR (L)# (5046664. 3-
* (4977793.7+(5017.8 *PSINEB)))**.5)72508.9
RKL(IP) = SATK(L)*(RMSTL(IP)/POR(L))**!!.76
CL(IP) = POR(L)/(5046664.3-
* (4977793.7+(5017.8*PSINE6>))**0.5
60TO 20
Soil type 3
203 CONTINUE
RMSTL (IP) = POR (D*. 89030+(POR (D* (27432013.6-
* (27015529.7 +(11765.B*PSINES)))»*.5)75882.9
RKL(IP) = SATK (D* (RMSTL (IP)/POR (L))**12. 80
CL(IP) = POR(L)/(27432013.6-
* (27015529.7 +(11765.8*PSINEB)))**0.5
SOTO 20
Soil type 4
204 CONTINUE
RMSTL (IP) = POR (D*. 88581+(POR (D* (333012588. 2-
* (327478085.9+(41202.4*PSINES)))**.5)/20601.2
RKL(IP) = SATK(L>*(RMSTL(IP)/PDR(L)>»*13. 60
CL(IP) = POR(L)/(333012588.2-
* (327478085.9+(41202.4*PSINEB)>)**0.5
SOTO 20
Soil type 5
205 CONTINUE
RMSTL (IP) = POR (D*. 88472+ (POR (L)*( 121120298. 2-
* (119063670.1+(24879.1*PSINES)))**.5)/12439.6
RKL(IP) = SATK (L>* (RMSTL (IP)/POR (L»**13. 78
CL(IP) « POR(L)/(121120298.2-
* (119063670.1+(24879.1*PSINE6)))**0.5
SOTO 20
Soil type 6
206 CONTINUE
RMSTL (IP) = POR (D*. 85500+(POR (D* (30352144. 1-
* (29479186.5 +(12887.2*PSINES)))**.5)/6443.6
RKL(IP) = SATK (D* (RMSTL (IP)/POR (L) )**17. 24
CL(IP) = POR(L)/(30352144.1-
* (29479186.5 +(12887.2*PSINEB)))**0.5
GOTO 20
Soil type 7
207 CONTINUE
RMSTL (IP) = POR (D*. 83734+(POR (D* (33601180.2-
* (32333107.1 +(13845.5*PSINEGM)**.5)/6922.7
RKL(IP) = SATK (D* (RMSTL (IP) /POR (L) )#*18. 50
CL(IP) = POR(L)/(33601180.2-
* (32333107.1 +(13845.5«PSINEB)))**0.5
SOTO 20
Soil type 8
208 CONTINUE
RMSTL (IP) = POR (D*. 80564+(POR (D* (69944841. 1-
* (65873969.3 +(20761.9*PSINEB)»**.5)/I0380.9
RKL(IP) = SATK (D* (RMSTL (IP)/POR (L))**20. 04
CL(IP) = POR(L)/(69944841.1
G-29
-------
* -(65873969.3 +(20761.9*PSINEG)>)**0.3
GOTO 20
C Soil type 9
209 CONTINUE
RMSTL (IP) = POR (D*. 54182+(POR (D* (347853. 3-
* (99106.0 +(2177.1 *PSINE6)))**.5)/1088.5
RKL(IP) = SATK (D* (RMSTL (IP)/POR (L»**23.80
CL(IP) = POR(L)/(347853.3-
* (99106.0+(2177.1 *PSINEG)))**0.5
GOTO 20
C Soil type 10
210 CONTINUE
RMSTL(IP) » POR(L)».54182+(POR(L)*(3567840.6-
* (1016504.7 +(6972.3 *PSINE6)))**.5)X3486.2
RKL(IP) = SATK(L)*(RMSTL(IP)/POR(L))**23.80
CL(IP) - POR(L)/(3567840.6-
* (1016504.7 +(6972.3 *PSINE6)))**0.5
GOTO 20
C Soil type 11
211 CONTINUE
RMSTL (IP) = POR (D* (-3.64) +(POR (L)*
* (1074124.4+(671246.7-(569.5*PSINEG)))»*.5)/284.7
RKL(IP) = SATK (D* (RMSTL (IP)/POR (L))*#25. 80
CL(IP) = POR(L)/(1074124.4+
* (671246.7 -(569.5 *PSINEG)))*#0.5
GOTO 20
C Soil type 12
212 CONTINUE
RMSTL(IP) * POR(L)
RKL(IP) = SATK(L)
CL(IP) = 0.0D0
END IF
20 CONTINUE
10 CONTINUE
C Compute geometric mean interblock conductivity
DO 120 L * 1,NUMEL
CDP CHANGE NEXT CARD FOR DOUBLE/SINGLE PRECISION
STARK(L) = SORT(RKL(L*2-1)*RKL(L*2)>
120 CONTINUE
C Compute mean RMQIST, RK, and C at nodes
DO 30 I - 2,NPM1
LI » 2*(I-1)
L2 = LI + 1
C(I) - (CL(L1)+CL(L2))/2.D0
RMOIST(I) = (RMSTL(LI)+RMSTL(L2))/2.D0
CDP CHANGE NEXT CARD FOR DOUBLE/SINGLE PRECISION
RK(I) = SQRT(RKL(L1)*RKL(L2>)
30 CONTINUE
C Use one node value at top and bottom
G-30
-------
C<1> = CL<1)
RMOIST(l) = RMSTL(l)
RK<1) = RKL(l)
C(NUMNP) = CL(NUMEL2>
RMOISTCNUMNP) - RMSTL(NUMEL2)
RK(NUMNP) = RKL
-------
APPENDIX H
SOILINER PACKAGE
H-l
-------
USING THE SOILINER PACKAGE
The SOILINER software package consists of three distinct modules:
(1) PRESOIL.BAS, the BASIC input preprocessor,
(2) SOILINER.EXE, the numerical model, and
(3) SGRAPH.WKS, the LOTUS 1-2-3 graphics macro.
The first two are linked together through the batch file SL.BAT. In addition,
the package includes 8 default data sets for a variety of clay liner
configurations, assorted utility files, and the FORTRAN source code for the
model. The entire package is stored on one double-density, double-sided,
5 1/4 inch floppy disk. The package was developed for use on the standard EPA
microcomputer configuration -- an IBM PC-AT equipped to display and print
LOTUS graphics. Equivalent software is available for use on a Compaq or an
IBM PC-XT, the only requirements being (a) a hard disk with LOTUS drivers
properly installed, and (b) graphics capabilities. A desirable option is a
math coprocessor (Intel 8087 for Compaq and XT, 80287 for AT), which
significantly accelerates execution time. This manual describes installation
and use of the SOILINER package, valid for all the above mentioned hardware
configurations.
***This manual references a number of keystrokes. If the user is asked
to type a command, only those letters within the quotes are to be typed.
When the user is asked to press a key (or sequence of keys) the text
within quotes usually represents a function key (i.e. "Return").
Simultaneous keystrokes are denoted by brackets (i.e.
-------
2.2 Type "SL" and press "Return". This invokes the batch file SL.BAT. At
this point, the user will be prompted for the first decision, whether to edit
input data sets or not. If editing is selected, the user must first choose
the flow scenario to be simulated. A second prompt will appear asking the
user whether the original, default data sets or the most recently revised data
sets are to be edited. Finally, the user is given a choice of editing the
control file, the soil properties file, and/or the initial conditions file.
Editing lasts until the user finally opts to exit the preprocessor, at which
time the SOILINER program is invoked.
If the edit option is not chosen, the user will be prompted to select the
default flow scenario to be simulated which, then invokes the SOILINER model.
If one chooses to run user-created input data sets instead, a final series of
prompts will appear requesting the appropriate data set names. Upon
successfully entering all four data set names, the simulation will proceed.
*** Note that the original default data sets have the extension .DBF and
cannot be changed, whereas the most recently edited data sets have the
extension .DAT.
Section 3 - EXECUTION OF THE MODEL
3.1 It is recommended that the user press upon exiting the
preprocessor, to print the screen activity of the simulation.
3.2 After exiting the preprocessor, the SOILINER program will execute. If
the steady state algorithm is chosen (LSTEDY = T in the control file), the
only screen display will be a brief summary of the results: (1) the number of
iterations, and (2) the error at the last iteration. If the number of
iterations exceeds the specified maximum number of iterations (ITER > MAXIT),
a warning will be issued and the user should check for convergence.
If the transient solution strategy is chosen (LSTEDY = F in the control file),
steady state information will be followed by time-stepping data (including the
particle's position) and the times to steady state and breakthrough.
Section 4 - USING SGRAPH
Due to the flexibility of the SOILINER graphics post-processor, SGRAPH,
this section of the manual is more detailed in scope.
Section 4.1 - Invoking SGRAPH
4.1.1 Get into the LOTUS directory by typing "CD\LOTUS" and pressing
"Return".
4.1.2 Enter LOTUS 1-2-3 by typing "LOTUS" and pressing "Return".
4.1.3 Select 1-2-3 from the menu by pressing "Return". At the prompt,
strike any key.
*** If you hear a beep and get an error message at this point, press
H-3
-------
"Return" and you should see a blank worksheet.
4.1.4 To retrieve the SOILINER worksheet, enter the 1-2-3 main menu by
pressing "/", then select "FILE", then "RETRIEVE", then "SOILINER" from the
consecutively displayed menus. The SOILINER worksheet (copied into the LOTUS
directory during installation) will, in turn, retrieve the SGRAPH worksheet
from C:\SOIL.
*** To make a selection from a menu, point to the desired selection using
the arrow keys, and press "Return". One can also select from a menu on
the very top line by typing the first letter of the desired option.
*** If SOILINER is not listed as an available worksheet, press "Esc",
then "/" to enter the 1-2-3 main menu, then select "FILE" and then
"DIRECTORY". Then type "C:\LOTUS" and press "Return", and perform step
4.1.4. again.
Section 4.2 - Creating Graphics
The following procedure will convert two sets of SOILINER data (pressure
and moisture) into graphics using SGRAPH. Following this procedure the first
time will provide familiarity with the SGRAPH COMMAND (CMD) MENU. Be sure
that the 1-2-3 and PRINTGRAPH device drivers are properly installed. Refer to
Section 4.3 when in trouble; Section 5 presents a breakdown of all SGRAPH menu
options.
4.2.1 Select "IMPORT" from the SGRAPH main menu.
4.2.2 Select "PRESSURE" from the next menu. A list of data files will
appear. Select the SOILINER pressure output file (PSI is the most recent)
which you wish to graphically display.
4.2.3 Sit back and wait. This process takes a few moments, and the screen
will jump and flash as the data is imported and prepared for graphical
presentation.
4.2.4 Select "MOISTURE" from the displayed menu, and a list of data files
will appear again. Select the SOILINER moisture output file (MOIST is the
most recent) which you wish to graphically display, and wait again.
4.2.5 At this point the worksheet contains two sets of data, the maximum
SGRAPH can handle at one time. Select "QUIT", then "GRAPH".
4.2.6 To graph the pressure data, select "PRESSURE".
4.2.7 Select "VIEW" to see the graph on the screen. Strike any key to
return to the options menu.
4.2.8 Assuming that the graph is satisfactory, it may be saved for later
printing. Refer to Section 5.0 to modify graphs using SGRAPH, or the LOTUS
1-2-3 manual to modify graphs using 1-2-3 commands directly. Select "SAVE",
or skip to Section 4.2.11 if no graphics hadcopy is desired.
H-4
-------
4.2.9 Enter a name for the graph being saved and press .
4.2.10 Select "QUIT" to return to the main graph menu.
4.2.11 To graph the moisture data, select "MOISTURE" and repeat instructions
4.2.7 through 4.2.10.
4.2.12 At this point both graphs are prepared for printing. Quit to READY
mode by selecting "QUIT" or pressing "Esc" three times (see section 4.3).
Press "/" for the 1-2-3 main menu, then select "QUIT", then "YES".
4.2.13 Now in the LOTUS Access System menu, select "PRINTGRAPH" to create
graphics hardcopies.
4.2.14 If PRINTGRAPH is not configured to read picture files from C:\SOIL,
refer to the 1-2-3 manual for instructions on specifying C:\SOIL as the
picture menu.
4.2.15 Select "SELECT", then the graphs to print. Making sure the printer is
on line, select "GO".
Section 4.3 - Error Correction Keys
Users experienced with 1-2-3 will be accustomed to hitting the "Esc" key
to "back out" of any unfamiliar territory they might run into in the 1-2-3
menu. This will work most of the time in SGRAPH, but not always. The
following keystroke sequences will get you out of almost any unpleasant
situations you might encounter in these steps - expect to become very familiar
with all of these as you use SGRAPH.
4.3.1 "Esc", "Esc", "Esc",...
First, press which stops SGRAPH. Then press "Esc" until
the word READY appears in the mode indicator in the top right corner of the
screen. Then press to get back into the SGRAPH menu. This is the
safest and most versatile escape mechanism, but has one quirk, in that SGRAPH
will sometimes hang on to the command. If the computer beeps
and you see "CTRL-BREAK" in the bottom left corner of the screen, try step
4.3.2.
*** If you should strike the "Break" (Scroll Lock) key without the "Ctrl"
key, the word "SCROLL" will appear in the bottom right of the screen. If
this happens, press "Scroll Lock" again to make this indicator disappear.
4.3.2 "Esc", "Esc",... "Esc"
Hit "Esc" until you are in familiar territory. If, in doing so, you find
yourself in READY mode, press . Or, if the mode indicator reads CMD
READY, press "Return", then . This sequence of key strokes should put
you back into the SGRAPH menu.
H-5
-------
4.3.3 QUIT (or Q)
QUIT is an option on all SGRAPH submenus. Select "QUIT" or type "Q" to
back up to a more general menu, or to get back to the READY mode. The QUIT
option is generally identical in function to the Esc key.
4.3.4
Should you select "QUIT" too many times, you lose the SGRAPH menu
altogether and wind up in READY mode. If you see the word "READY" in the top
right hand cell of the screen, press to reenter the SGRAPH menu
system.
Section 5 - THE SGRAPH MENU SYSTEM
The SGRAPH menu system is outlined in the table below. Each row of this
table represents a single selection/command sequence. All commands are marked
by asterisks and described in Sections 5.2 through 5.13. The remaining
selections will route the user to a new menu.
5.1 The SGRAPH Menu
SELECTIONS AND COMMANDS
LEVEL I LEVEL II LEVEL III LEVEL IV LEVEL V
IMPORT * PRESSURE
IMPORT * MOISTURE
IMPORT * CLEAR
IMPORT QUIT
GRAPH
GRAPH
GRAPH OPTIONS XYSCALES XSCALE
GRAPH OPTIONS XYSCALES XSCALE
GRAPH OPTIONS XYSCALES XSCALE
GRAPH OPTIONS XYSCALES YSCALE
GRAPH OPTIONS XYSCALES YSCALE
GRAPH OPTIONS XYSCALES YSCALE
GRAPH OPTIONS XYSCALES QUIT
GRAPH
GRAPH
GRAPH
GRAPH
GRAPH
QUIT
All QUIT selections lead back to the previous menu, except for the Level I
QUIT selection, which leads to the 1-2-3 ready mode. To enter the Level I
menu from READY mode, press .
* PRESSURE
* MOISTURE
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
OPTIONS
QUIT
XYSCALES
XYSCALES
XYSCALES
XYSCALES
XYSCALES
XYSCALES
XYSCALES
* VIEW
* NAME
* SAVE
QUIT
* AUTOMATIC
* MANUAL
QUIT
* AUTOMATIC
* MANUAL
QUIT
H-6
-------
5.2 PRESSURE (Import Menu)
This command takes pressure data from a SOILINER output file and puts it
into a predesignated range on the worksheet. Using this command will
automatically delete any pressure data currently on the worksheet, but must be
executed before the pressure data can be graphed.
5.3 MOISTURE (Import Menu)
Identical to 5.2, except that it is used for importing moisture data.
5.4 CLEAR (Import Menu)
This command erases all pressure and moisture data from the worksheet,
leaving a clean sheet. Use this command for aesthetic purposes only - it is
not necessary to execute this command before importing additional sets of
pressure or moisture data.
5.5 PRESSURE (Graph Menu)
This command takes previously imported pressure data and puts it into a
basic graphic format, ready for viewing. Note that while SGRAPH accomodates
both pressure and moisture data simultaneously, it normally will store only
one graph at a time. This function converts the pressure data into the
current graph. Exercise all desired graphics options (XYSCALES, VIEW, NAME,
SAVE) after using this function and before graphing a new data set.
5.6 MOISTURE (Graph Menu)
Identical to 5.5 but for moisture data.
5.7 XYSCALES XSCALE Automatic (Options Menu)
This step tells SGRAPH to automatically set the X (horizontal) scale on
the current graph. Note that 1-2-3 will normally do this by default. Use
this step only if you have tried adjusting the X scale manually and prefer the
original format.
5.8 XYSCALES XSCALES Manual (Options Menu)
This step allows you to set the X scales manually. It will request a
lower and an upper limit - press RETURN after entering each. Remember that
you can always use step 5.7 if this function leads you astray.
5.9 XYSCALES YSCALE Automatic (Options Menu)
Same as 5.7, except that it sets the Y (vertical) scale.
H-7
-------
•;.]() XYSCALES YSCAhE Manual (Options Menu)
Same as 5.8, except that it enables you to set the Y (vertical) scale.
5.11 VIEW (Options Menu)
Use this function to view the current graph on the screen. Press any key
to get back to the OPTIONS menu.
5.12 NAME (Options Menu)
This is one of the links between SGRAPH and the 1-2-3 main menu. Use it
if you are somewhat familiar with 1-2-3 graphics and want to store more graphs
on the worksheet than just the current one. Enter a new name each time you
use this function, since SGRAPH contains no provision to avoid overwriting
graphs already named on the worksheet.
5.13 SAVE (Options Menu)
This is the major link between SGRAPH and Lotus's PRINTGRAPH. It saves
the current graph in the directory C:\SOIL for future printing using
PRINTGRAPH. See the 1-2-3 manual for a detailed explanation of PRINTGRAPH.
*** Remember that SAVE commands in the SGRAPH menu do not contain any
provisions to avoid overwriting of old files. Type a new name for each
new file you create.
H-8
-------
DEFAULT DATA SETS
Included in the SOILINER package are default data sets for eight flow
scenarios. Associated with each scenario are the four required data sets, as
listed below:
1. Clay liner only (60 cm)
CTRC60.DEF
GRDC60.DEF
PRPC60.DEF
INTC60.DEF
3. Clay liner only (180 cm)
CTRC180.DEF
GRDC180.DEF
PRPC180.DEF
INTC180.DEF
Clay liner (60 cm) with
underlying native soil
(300 cm)
2. Clay liner only (90 cm)
CTRC90.DEF
GRDC90.DEF
PRPC90.DEF
INTC90.DEF
4. Clay liner (60 cm), drainage
layer (30 cm), and clay
liner (90 cm)
CTRCSC.DEF
GRDCSC.DEF
PRPCSC.DEF
INTCSC.DEF
6. Clay liner (90 cm) with
underlying native soil
(300 cm)
CTRCS60.DEF
GRDCS60.DEF
PRPCS60.DEF
INTCS60.DEF
7. Clay liner (180 cm) with
underlying native soil
(300 cm)
CTRCS180.DEF
GRDCS180.DEF
PRPCS180.DEF
INTCS180.DEF
CTRCS90.DEF
GRDCS90.DEF
PRPCS90.DEF
INTCS90.DEF
8. Clay liner (60 cm), drainage
layer (30 cm), clay liner
(90cm), with underlying
native soil (300 cm)
CTRCSCS.DEF
GRDCSCS.DEF
PRPCSCS.DEF
INTCSCS.DEF
Note that for each set of input data files, CTR, GRD, PRP, and INT
represent the control, grid, soil properties and initial conditions files
respectively.
H-9
------- |