c/EPA
United States
Environmental Protection
Agency
EPA/530-SW-84-001
APRIL 1984
c,
Solid Waste
Procedures for
Modeling Flow
Through Clay
Liners to Determine
Required Liner
Thickness
Draft
Technical Resource
Document For Public
Comment
-------
PROCEDURES FOR MODELING FLOW
THROUGH CLAY LINERS TO
DETERMINE REQUIRED LINER THICKNESS
EPA/530-SW-84-001
This report was prepared for the
Office of Solid Waste
under Contract No. 68-02-3168
U £ r- -.,.-. •• '• - . Agency
r:
2 J.!:t': L, ii UL:"C...-. „„•;**"'
Chicago, Hiinoi5 60604
U.S. ENVIRONMENTAL PROTECTION AGENCY
1984
-------
DISCLAIMER
This report was prepared by Daniel J. Goode and Patricia A. Smith of
GCA Corporation, Technology Division, Bedford, Massachusetts under
Contract No. 68-02-3168. 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. There-
fore it does not necessarily reflect the views or policies of the Agency.
Mention of trade names or commercial products does not constitute endorse-
ment or recommendation for use.
-------
ACKNOWLEDGMENTS
The authors, Daniel J. Goode and Patricia A. Smith, express their
appreciation to Dr. P. Christopher D. Milly who provided guidance and
review during this project and who computed verification analytical solu-
tions for model testing. The burden of producing this technical document
in a readable form was carried by GCA's fine publications staff, to whom
the authors are indebted.
11
-------
FOREWORD
The Environmental 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.
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 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 single soil liners. It includes a perfor-
mance simulation model that is based on numerical techniques
recommended in guidance.
John H. Skinner
Director, Office of Solid Waste
U.S Environmental Protecton Agency
ill
-------
CONTENTS
Figures v
Tables vii
1. Introduction 1
2. Conceptual Model of Flow through Liner 2
2.1 Description of Physical Problem 2
2.2 Mathematical Statements 2
2.2.1 Governing Equations of Vertical Unsaturated Flow ... 2
2.2.2 Boundary Conditions 4
2.3 Soil Moisture Characteristics 6
2.4 Unsaturated Hydraulic Conductivity 12
2.5 Verification of Conceptual Model 14
3. Determination of Soil Hydraulic Properties 20
3.1 Introduction 20
3.2 The Soil Moisture Characteristic 20
3.3 Unsaturated Hydraulic Conductivity 21
3.3.1 Laboratory Methods 21
3.3.2 Field Methods 22
4. Numerical Simulation of Unsaturated Flow 23
4.1 Finite Difference Method 23
4.1.1 FDM Spatial Difference Approximations 23
4.1.2 FDM Temporal Derivative Approximations 25
4.1.3 FDM Application to Vertical Unsaturated Flow 28
4.2 Finite Element Method 28
4.2.1 FEM Spatial Approximation 28
4.2.2 FEM Applications to Vertical Unsaturated Flow 31
4.3 Discretization of Flow Domain 31
4.3.1 Grid Design 31
4.3.2 Time Stepping 32
4.4 Soil Properties in Numerical Models 32
4.4.1 Tabular Interpolation 32
4.4.2 Functional Relationships 37
4.5 Example Compute Program 39
4.5.1 Numerical Technique 39
4.5.2 Soil Property Representations 42
4.5.3 Verification 42
5. Soil Liner Design 50
5.1 Introduction 50
5.2 Discretization of the Flow Domain 50
5.3 Soil Properties from Data and Models 50
iv
-------
CONTENTS (continued)
5.4 Simulation 52
5.4.1 Infiltration 52
5.4.2 Breakthrough 52
5.5 Adjusting Liner Specifications 56
5.6 Summary 56
6. Summary 59
7. References 60
Appendices
A Review of the Transit Time Equation for Estimating Storage Impoun-
ment Bottom Liner Thickness A-l
B Partial List of Available Unsaturated Flow Models B-l
C Listing of Example Computer Program Soiliner C-l
D Example Input and Output for Soiliner Model D-l
E Computer Program for Gardner's Analytical Solution E-l
-------
FIGURES
2.1 Flow domain for liner breakthrough 3
2.2 Schematic of initial capillary liquid pressure distribution ... 5
2.3 The effect of texture on soil moisture characteristic 7
2.4 A piecewise linear relation between moisture content and the
logarithm of suction 10
2.5 Sandy and clayey soil moisture characteristics computed from a
continuous functional relationship 11
2.6 Schematic of unsaturated hydraulic conductivity for a sand and a
clay soil 13
2.7 Experimental and theoretical moisture content profiles for
Columbia silt loam 15
2.8 Computed (FDM) and measured field soil moisture profiles before
and after 9 1/3 hours infiltration 16
2.9 Comparison of compated infiltration into Yolo light day with
ponding by FEM and quasi-analytic methods 17
2.10 Computed dimensionless pressure head as a function of depth and
time for infiltration from an axisymraetric infiltrometer. ... 19
4.1 Finite difference method spatial discretization grids 24
4.2 Finite difference method temporal discretization for node i at
time level n + 1 showing explicit and implicit relationships. . 26
4.3 Finite element method discretization and linear interpolation
functions 29
4.4 Two vertical grids with variable subdomain sizes 33
4.5 Interpolation of specific moisture capacity from table of values
0.5 pF apart 35
4.6 Interpolation of moisture content from table of values
O.5 pF apart 36
4.7 Interpolation of unsaturated hydraulic conductivity from values
0.5 pF apart 38
4.8 SOILINER solution procedure flow chart 43
4.9 Comparison of solutions for steady state vertical flow upward
from a water table 46
4.10 Comparison of Philip's quasi-analytic solution and SOILINER
regular grid solution for infiltration into Yolo light clay
under ponding 47
4.11 Comparison of Philip's quasi-analytic solution and SOILINER
graded grid solution for infiltration into Yolo light clay
under ponding 48
vl
-------
FIGURES (continued)
5.1 Soil properties of hypothetical clay and sand for liner design. . 51
5.2 Simulated pressure versus time at several points in liner with
initial moisture content of about 0.25 54
5.3 Simulated moisture content profile in liner and site soil with
initial moisture content of about 0.25 55
5.4 Simulated pressure versus time at several points in liner with
initial moisture content of about 0.30 57
5.5 Simulated moisture content profiles in a three-layer liner and
underlying site soil 58
vii
-------
TABLES
5.1 Comparison of Simplified Models for 5-year Liner.
Page
53
viii
-------
1. INTRODUCTION
The Part 264 Subpart K regulation and associated guidance (47 FR 32274)
broadens the legal storage options available to owners and operators of
hazardous waste facilities insofar as allowance is made for the use of
adequate single soil liners to create surface impoundments for temporary
storage of wastes. Upon closure of the impoundment, all wastes and all parts
of the liner permeated by wastes are to be removed for treatment, storage, or
disposal at a RCRA hazardous waste management facility.
In order to assist the design and engineering of adequate single clay
liners, GCA prepared this report which details the procedures for determining
adequate thicknesses of single soil liners. An example is presented of a
performance simulation model which is based on numerical techniques
recommended in draft guidance.
This report includes a conceptual model of vertical unsaturated flow
through porous media and the governing differential equations (Section 2);
referenced methods for conducting field and laboratory tests to generate the
soil property data required for model use and liner design (Section 3); and an
example of computer-assisted procedures for modelling flow through natural
clay liners. Finite difference and finite element techniques were selected
for analysis of the models nonlinear governing equations (Section 4). These
numerical procedures have sufficient flexibility to be able to incorporate
many different soil types, soils with spatially varying properties, and
temporal variations. The utility of the model in design evaluation and the
application of a GCA-developed finite difference model are demonstrated
through simulation of a hypothetical liner's performance.
-------
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
natural soil liner of thickness d (L); a layer of underlying site soil, which
may or may not be saturated; 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 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 reaching ground water during the impoundment's design life.
2.2 Mathematical Statements
2.2.1 Governing Equation of Vertical Unsaturated Flow—
The governing equation for one-dimensional unsaturated flow in the
vertical direction can be written [see Bear, 1979, p. 214]:
, - - * - - , (2-1)
3z L 3Z J
in which ty&$ -z [L] is matric potential or capillary pressure head, where
<)> [L] is piezometric head; 0 [-] is volumetric moisture content; K(i^) =
I^djj) Ks [LT~M is unsaturated hydraulic conductivity, where Kr(^) [-]
is relative hydraulic conductivity, and Ks [LT~M is saturated hydraulic
conductivity; z [L] is the vertical coordinate, positive upward; and t [T] is
time. This equation is developed by consideration of conservation of mass.
The first term represents the change in storage of liquid mass and the second
term is mass flux divergence, or the change in flux over space. The second
term contains fluxes due to the matric potential gradient (K(i^) 9\j;/Sz) and
gravitational potential (K(^) 9z/9z). The flux term is developed from the
generalization of Darcy's law for water flow in porous media,
(2-2)
-------
DATUM
IMPOUNDED
LIQUID
/
/
SOIL
LINER
MATERIAL
-z,
\/
NATURAL
FILL
'"SATURATED,*.
o;/. ZONE • .V*
• •„»•••*;
n ° « O »*^ <
Figure 2-1. Flow domain for liner breakthrough.
-------
in which q [LT 1] 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-1), derived for unsaturated flow, is also applicable to
modeling temporarily or permanently saturated soil zones. In that case, ^ is
known as the pressure head, d 6/d becomes the specific storativity, and K.OJJ)
is equal to Kg.
2.2.2 Boundary Conditions —
At the top of the liner, z=0, the liquid pressure 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:
$ = $ « ~ z at z=0 (2-3)
z=0
Since the piezometric head is constant in space (hydrostatic) in the ponded
liquid,
* z=0 = *z=h = h£ (2~4)
and, since z=0, (2-3) becomes:
^=(l)z=0 -0 = h£ at z=0 (2-5)
This is the fixed pressure Dirichlet boundary condition applied at the top of
the soil column.
By definition, the matric potential at the water table is equal to zero;
thus at the column bottom, z = zw, the Dirichlet boundary condition is:
4> = zw or 4* = 0 at z = z (2-6)
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]:
4* ='IJ at z = -d (2-7)
s
where tp and 4igj refer to liner and site soil respectively. This is simply a
condition of pressure continuity. This condition is incorporated into the
governing equation (2-1) and is implicitly satisfied.
-------
t-
o
00
z
o
LU
_l
LU
IMPOUNDMENT
HYDROSTATIC —
PRESSURE
SOIL
LINER
V
CONSTANT
INITIAL
PRESSURE
UNSATURATED
NATURAL
FILL
I
HYDROSTATIC
PRESSURE
I
I
SATURATED
ZONE
I
T
NEGATIVE (SUCTION) (POSITIVE
LIQUID PRESSURE
Figure 2-2. Schematic of initial capillary liquid pressure distribution.
-------
2.3 Soil Moisture Characteristics
Liquid is held in the pore space of an unsaturated porous medium by
capillary and adsorptive forces [Hillel, 1971, p. 57]. In a soil column at
static equilibrium, these surface tension and adsorptive forces among the
liquid, the solids, and the air exactly oppose the force of gravity pulling
the liquid down. In general, the smaller the pore size, the more weight these
forces can support. For the liquid to be at static (no flow) equilibrium, the
piezometric head must be constant ( ^ = ^ +z = constant) over the entire
column. Thus, the matric potential decreases linearly with height:
*=<(> - z (2-8)
constant
At and below the water table = zw and above the water table
* = z - z (2-9)
w
This static pressure profile is shown in Figure 2-2.
Figure 2-3 shows typical soil moisture characteristic curves for a sandy
and a clayey soil. Clays hold more liquid at lower pressures (or, in a
hydrostatic situation, higher elevations above the water table) because of
their smaller pore size. Although the porosity of a clay may be higher than
that of a sandy soil, the individual pores are much smaller.
To evaluate the storage term of the governing equation (2-1), we must
define the soil specific moisture capacity C( ijj) [L~l] :
which is the derivative of the moisture retention function
0 =6 (*) (2-11)
that describes the relationship between moisture content and matric potential,
shown in Figure 2-3.
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:
n
(n -
JiL
*b
*>*b
* < \ (2-12)
— b
in which or [-] is residual (nonreducible) moisture content, ^ [L] is
bubbling or air entry pressure at which air first enters the draining column,
n [-] is porosity, and A is a fitted parameter. The experimental data to
develop this model included pressures of 0 _> i|> 2. "500 cm. The corresponding
specific moisture capacity is:
-------
•s-
I
o>
o
UJ
o
Q.
O
or
-CLAYEY SOIL
SANDY SOIL-
Q- WATER CONTENT
Figure 2-3. The effect of texture on soil moisture characteristics
(after Hillel, 1971).
-------
c(40 =
(2-13)
This function is discontinuous at saturation, 41 = 4^. Clapp and Hornberger
[1978] amended this model to include gradual air entry near saturation and to
provide a continuous relation.
King [1965] fit a hyperbolic function to moisture retention data for
several soil types with pressures down to 41 = -100 cm:
= n<5
cosh [(4V4- ) +e ] -y
o
cosh
+E ]
(2-14)
in which 5, ^Q, [J, c, and y are fitted parameters. This model was used by
Gillham et al. [1976] to provide data for numerical computations.
Rogowski [1971] proposed a simple model requiring few input parameters.
This model is:
n + a In (4> -u;+l)
b
(2-15)
= n
in which
= (9, - n) / In ( ^ - ty
ID b 15
1)
(2-16)
where 4^5 = 1.5 x 10^ cm (15 bars) and 6^5 is moisture content when
fy = 4)15- Again, this model is discontinuous at saturation ( 4) = 4*^) but
only requires three measured parameters: n, 4^, and 0^5. This model
performed well for several soil types, including a clay.
McQueen and Miller [1974] proposed a linear relationship between pF =
log (-41). with 41 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:
6(pF)
- In [exp |M (a -s pF)> + exp < M
PF)i]
- ~ In [exp
pF)|
(2-17)
exp
-------
in which a , a , s , s and 9 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:
9(40 = a2 - s2 log (-40 (2-18)
and the specific moisture capacity is approximately
£n(lO)
(2-19)
Since (2-17) has continuous slope, the exact form of (2-19), for all ty, is
also continuous and can be written:
-1
(2-20)
-1
_
M1
where:
M a,
Pl =
q1 = - M s- log (e)
M a.
P2 -e
(2-21)
M'a,
P3 -e
q = -M's log (e)
M1
P4 =e
Equations (2-17) and (2-20) are plotted for sandy and clayey soils in
Figure 2-5.
The moisture content of a soil is not a unique function of matric
potential; this relationship varies widely between wetting and drying
processes. This effect is called hysteresis and is important in determining
redistribution of moisture in a soil column when some spaces are being wetted
and later dried [Hillel, 1971]. During infiltration under ponding, the effect
of hysteresis is small if the wetting history of different points in the soil
is similar and all pores are wetting [see Morel-Seytoux, 1973, p. 176; Dane
and Wierenga, 1975].
-------
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).
10
-------
a.
15 .20 .25 .30 .35
MOISTURE CONTENT
b.
(U.
I
n.
i
a
_J U>-J
(a —
i
YOLO
SflMD
.0
Figure 2-5.
1.0
1.5
PF
S.Q
2.5
3.0
Sandy and clayey soil moisture characteristics computed from
a continuous functional relationship developed by Milly and
Eagleson (1980), a. moisture content and b. log of specific
moisture capacity (cm~l) vs. PF = log(-^), ^in cm.
11
-------
2.4 Unsaturated Hydraulic Conductivity
Darcy's law for flow in porous media states that the flux rate is
rtional to the total piezometric head gradient times the hydraulic
proportional
conductivity [Bear, 1979]
q = -Kty) (2-22)
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. This assumption can be
relaxed in a more general approach to the problem.
For a given soil, the hydraulic conductivity is dependent on the moisture
content or matric potential in the soil. As a function of moisture content,
6, the hydraulic conductivity function shows little hysteresis [Bear, 1979;
Mualem, 1976] and since, as statpH above, Q(^) hysteresis for this
infiltration problem is ignored, the hydraulic conductivity as a function of
matric potential can also be determined uniquely. Figure 2-6 shows a
schematic of typical relationships between hydraulic conductivity and matric
potential for sandy and clayey soil.
Numerous functional relationships have been proposed for the unsaturated
hydraulic conductivity. These models include
Gardner [1958]: K(); (2-24)
c
s
—m
Brooks and Corey [1964]: K( d9 (2-28)
. w
12
-------
K
sat
Ksat
X.
o
Figure 2-6. Schematic of unsaturated hydraulic conductivity for a sand
and a clay soil.
13
-------
in which yw is the specific weight of water. These formulas (2-26 - 2-28)
require use of cgs units. Equation (2-28) "represents the amount of work
required to drain a unit volume of a saturated soil" [Mualem, 1978]. Thus,
(2-26) is dependent on the moisture characteristic curve. This model was
shown to improve the prediction of unsaturated hydraulic conductivity for
fine-grained soil [Mualem, 1978] relative to previous techniques.
There are several integral techniques which, as (2-26 - 2-28) above,
derive the hydraulic conductivity function from the moisture characteristic
curve [see Mualem 1976; Elzeftamy 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 by a particular technique.
2.5 Verification of Conceptual Model
Several investigators have shown that the presented governing equation
for vertical unsaturated flow (2-1) accurately described the infiltration of
moisture into soil. This verification has included comparison of laboratory
and field experiments to analytical, finite difference, finite element, and
other numerical solutions, all based on the conceptual model (2-1).
The agreement between Philip's [1958; 1969] quasi-analytical solution and
experimental water content profiles reported by Davidson et al. [1969] is
shown in Figure 2-7 [after Swartzendruber, 1969]. Youngs [1957] reported
good agreement between theoretical and observed profile shape in laboratory
experiments. Nielsen et al. [1961] compared solutions of (2-1) to the results
of field studies and showed fair agreement.
The governing equation for vertical infiltration (2-1) has been solved
numerically by many authors. Green et al. [1970] showed good agreement
between a finite difference model of (2-1) and field measurements (see
Figure 2-8). Giesel et al. [1973], Elzeftawy and Dempsey [1976], Haverkamp
et al. [1977], and Ragab et al. [1982] compared results of laboratory
experiments and finite difference solutions. Kunze and Nielsen [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. The finite element technique was applied to the horizontal
unsaturated flow equation with functional soil properties and compared to
laboratory measurements by Hamilton et al. [1981]. Milly [1982] and Segol
[1976] showed the agreement between results from finite element models and
Philip's [1958; 1969] quasi-analytical infiltration solution (Figure 2-9).
As described above, the analysis herein considers only vertical flow.
Practically speaking, the liner thickness at a typical impoundment will be two
to three orders of magnitude smaller than the horizontal dimensions.
Horizontal flow may be significant near the edges of the impoundment, but this
effect should slow vertical movement. Jeppson et al. [1975] compared a
numerical solution of the one-dimensional governing equation (2-1) to a
14
-------
-30
"--40-
CALCULATED
FOR 1-07 hr.
EXPERIMENTAL
POINTS FOR'-
A - 1.07 hr.
D - 3.77 hr.
O - 7.78hr.
CALCULATED
FOR 3.77 hr.
CALCULATED
FOR 7.78hr.
-70-
0.1 0.2 0.3 0.4 0.5
MOISTURE CONTENT, 9
Figure 2-7. Experimental and theoretical moisture content profiles for
Columbia silt loam (after Swartzendruber,(1969), data of
Davidson et al. (1969)).
15
-------
UJ
u.
oc
o
z
<
UJ
03
Q.
UJ
0
10 20 30 40
PERCENT MOISTURE BY VOLUME
Figure 2-8. Computed (FDM) and measured field soil moisture profiles before
and after 91/3 hours infiltration, constant properties (after
Green et al., 1970).
16
-------
§b.20
MOISTURE CONTENT
0.30
•
0.40
I
2xlQ5 SECONDS
PHILIP (1958, 1969) QUASI-ANALYTIC
MILLY (1982) FEM
Figure 2-9. Comparison of computed infiltration into Yolo light clay with
ponding by FEM and quasi-analytic methods (after Milly, 1982).
17
-------
numerical solution of two-dimensional axisymmetric infiltration from a finite
diameter circular area. As shown in Figure 2-10, the wetting front at the
centerline of the axisymmetric solution advances slightly slower than the
wetting front from the one-dimensional solution. Thus, in homogeneous soils,
the one-dimensional analysis in conservative. For anisotropic layered soils,
with horizontal hydraulic conductivities higher than vertical, Siegel and
Stephens [1980], showed that lateral spreading can greatly reduce the vertical
rate of moisture front movement, and the one-dimensional analysis is even more
conservative.
18
-------
CENTER LINE AXISSYMETRIC
ONE-DIMENSIONAL VERTICAL
-5 -10 -15
PRESSURE HEAD, ^=
-20
-25
Figure 2-10.
Computed dimensionless pressure head as a function of depth
and time for infiltration from an axisymetric infiltrometer
(after Jeppson et al., 1975).
19
-------
3. DETERMINATION OF SOIL HYDRAULIC PROPERTIES
3.1 Introduction
The design of an adequate impoundment liner is grounded not only on the
conceptual model of flow developed but also on accurate assessment of the soil
properties of the liner and of the underlying site soil. Two such properties
are particularly relevant to the model's execution. These are the moisture
characteristic and the unsaturated hydraulic conductivity. This section
defines these properties and discusses their empirical determination.
Although well-documented and cost-effective methods exist for laboratory
determination of soil hydraulic properties, the accuracy of laboratory-obtained
values should be scrutinized (Roberts, 1982; Olson and Daniel, 1981). Of
foremost concern is the use of soil or permeant samples which may not be
representative of field conditions. Sampling and handling techniques may
significantly alter soil microstructure. In addition, field investigation and
statistical sampling, if inadequate, may fail to document macroscopic features
(i.e., fissures, root holes, sand seams or lenses) that significantly
influence local hydraulic properties. Finally, soil tests should be conducted
using a permeant similar to that to be impounded in order to include effects
of chemical and physical interactions in the values assessed for hydraulic
properties. Any additional cost associated with the adoption of in situ
measurement techniques should be weighed against the cost of using inaccurate
empirical parameters in the liner design. This section will refer to and
compare both field and laboratory testing procedures.
3.2 The Soil Moisture Characteristic
The soil moisture characteristic, described in Section 2.3, expresses the
relationship between water content of a soil and negative pressur2 or
suction. As sufficient negative pressure is applied—that is, as an outward
force is exerted and exceeds the capillary force which retains water in soil
pores—these pores begin to drain. Increasing suction causes drying as
smaller and smaller interstices empty. The moisture characteristic, a
hysteretic curve, illustrates that water is harder to withdraw from soil once
wet than it is to absorb once dry.
Hillel (1980) discusses methods that can be used in the laboratory or
in situ to measure soil moisture content continuously or intermittently.
These include neutron scattering, gamma-ray absorption, and electrical
resistance of porous blocks; the former two are suitable for use in field
studies. Once equipment is installed and calibrated, testing by these methods
requires considerably less time and effort than laboratory oven-drying
techniques which may yield erroneous results.
Tensiometers, pressure plates, and thermocouple psychrometers, each
accurate over different ranges of suction, can be used to monitor suction
continuously. Tensiometers are generally limited to suctions less than
0.9 atm in order to avoid cavitation of air bubbles in the device. The
20
-------
pressure plate is a modification of the tensiometer design capable of
accurately measuring suctions from 1 to 15 atm. While inaccurate for suctions
lower than 1-2 atm, psychrometers, by measuring the relative humidity of pore
air, can be used to determine suctions up to about 80 atm. When one of these
methods are used to obtain readings over time in conjunction with moisture
content readings, the wetting moisture characteristic for a design liner and
waste liquid can be constructed.
3.3 Unsaturated Hydraulic Conductivity
Hydraulic conductivity, K(ijj), described in Section 2.4, is a highly
variable soil property exhibiting a range of values over 14 orders of
magnitude from gravel to clays. The presence of organic solvents has been
documented as increasing the conductivity of clays (Haxo, 1976), sometimes
with no apparent tendency to reach maximum values (Brown & Anderson, 1982).
Although many techniques have been developed to measure hydraulic conduc-
tivities over this range, only a few are feasible for measurement of
conductivity in the type of unsaturated fine-grained soils that are suitable
for liners.
3.3.1 Laboratory Methods—
Olson and Daniel (1981) selected two tests for determination of
conductivity in unsaturated fine-grained soils.
A. The instantaneous profile test (Weeks and Richards, 1967, Watson,
1966). A soil sample is arranged in a column at equilibrium. By
imposing a constant or time-dependent suction or fluid flow,
unsteady-state seepage develops within the sample. Probes inserted
at different depths in the column measure water content and suction
values over time. A steady-rate of inflow test requires about
2 weeks to run.
B. The pressure plate method (Klute, 1965). This method entails
placing a soil sample on a saturated porous plate in a pressure
vessel. After a known air pressure is applied, water pressure in
the plate is kept at atmospheric pressure and the system allowed to
equilibrate. Then, by steps, the air pressure in the vessel is
abruptly altered, which causes excessive pore water pressure and the
drainage of water through the porous plate.
Olson and Daniel (1981) also discuss sources of error in laboratory
conductivity tests of partially saturated soils. They recommend as the best
laboratory techniques:
(1) At 0-0.9 atm of suction, the instantaneous profile method with
tensio.ne.ers;
(2) At 2-80 atm of suction, the instantaneous profile method with
psychrometers;
(3) At 1-15 atm of suction, pressure plate outflow.
21
-------
3.3.2 Field Methods--
Olson and Daniel (1981) suggest that although in situ unsaturated soil
field testing methods have only recently been developed, there are two
techniques available that can be used for the case of ponded fluids. Roberts
(1982) also presents an overview of these techniques, and related in situ
methods are described in Method 9100 of EPA guidance.
A. The instantaneous profile method has been adapted for field use. In
a diked test plot, probes for measuring water and suction pressure
are inserted at selected depths. The plot is flooded with the test
permeant and, after the permeant has soaked into the ground, covered
to prevent evaporation. If an appropriate moisture characteristic
has been previously constructed, moisture content probes or suction
probes alone can be used and the other value inferred from the
wetting branch of the characteristic (cf. Baker e_t_ aj^. , 1974; Rose
and Krishman, 1967). Testing times may be long in relatively
impervious soils.
B. Infiltration through impeding layer (Gardner, 1970; Hillel and
Gardner, 1970; Bouma, e_t^ aj^. , 1971). A column of soil is isolated
by pushing a thin-walled tube into the soil and then capped with a
relatively impervious porous stone or membrane. (This procedure
thus also may be conducted in the laboratory). Water is ponded on
top of the stone and a small constant head maintained to develop
steady-state seepage. Rate of flow through the stone is measured
with a Mariotte bottle. Tensiometers inserted into the soil are
used to measure suction and to confirm a hydraulic gradient of unity
directly beneath the stone. Given this unit gradient the hydraulic
conductivity is equivalent to the infiltration rate. Olson and
Daniel (1981) point out that for fine-grained soils, not only may it
be impractical to wait for steady state seepage to develop but the
stone impedance would have to be so great that accurate flow rate
measurements would be difficult to obtain.
22
-------
4. 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 nonhomogenicity of soil types
and 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. Available
computer-programs range from those treating basic flow [Bruch, 1975, Johnson,
et al., 1982] to coupled heat and moisture transport with vapor transport and
moisture hysteresis [Milly, 1982]. Using such programs, it is possible to
predict the performance of specific liner designs, with parameters from field
and laboratory tests.
4.1 Finite Difference Method
Finite Difference Methods (FDM) have been successfully applied to a wide
range of problems in groundwater flow and soil physics. 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 terms 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 terms. 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 up into a row of short vertical segments. This row of
segments is called a grid. For mesh-centered grids, the values of matric
potential (^) are evaluated at nodes which are located at the ends of each
segment (Figure 4-la). For a block-centered grid, these nodes are located at
the center of each segment (see Figure 4-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.
4.1.1 FDM Spatial Difference Approximations—
The spatial domain of the vertical unsaturated flow equation (2-1) 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. This domain is
discretized as shown in Figure 4-1 in which subscript i represents the node
number. Values of the state variable (matric potential, ^), and soil
properties (hydraulic conductivity, K(^), and moisture capacity, C(
-------
o. MESH CENTERED
b. BLOCK CENTERED
TOP
OF DOMAIN'
SOIL 1
SOIL
BOUNDARY
SOIL2
BOTTOM
OF
DOMAIN
NODE 1
NODE 2
NODE
NUMEL
NODE
NUMEL-H
SUBDOMAIN 1
SUBDOMAIN 2
SUBDOMAIN
NUMEL
NODE 1
2
NODE NUMEL
NUMEL = number of subdomoins
or elements
Figure 4-1. Finite difference method spatial discretization grids.
24
-------
G£ are functions of ^j_. These values are used to approximate the
derivatives in (2-1). The first term in (2-1), representing pressure driven
flux, can be replaced by:
"JN~~~ **• T~
3z 3z
v
i-L/2
1_
Az
in which Az =
Az i-1/2 Az
is the constant node spacing and
are the arithmetic average hydraulic conductivities between nodes.
gravitational flux term can be expressed as:
(4-1)
(4-2)
The
3KCIO
9z
K.
i+l
- K.
1-1
(4-3)
2Az
Equations (4-1) and (4-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. With constant node
spacing (4-3) is second order accurate. This means that if conductivity K(^)
is a liner function of z and z^ only, then (4-3) is exact.
4.1.2 FDM Temporal Derivative Approximations—
The temporal domain of the vertical unsaturated flow equation (2-1) is
time, t > 0, after infiltration into the liner begins. Time is broken down
into time level subdomains as shown in Figure 4-2 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: tyn £, Kn £, and Cn ^. The storage term in (2-1) can
be written:
r
C
n+u
It
n+1
(l-oi) |C
n+cx
At
(4-4)
in which Cn+a = acn+l + (1-a) Cn; 0 _^ a _< 1 is a temporal weighting
parameter, and At = tn+l - tn is the time step size. Subscripts have been
dropped for convenience. If a = 0, (4-4) is
n
2ft
(4-5)
which is forward differencing. When solving for \j;n+l, all other terms of
(4-5) are known from the last time step and thus 4in+l is solved explicitly.
If X=l, (4-4) is:
n+1
At
(4-6)
25
-------
TIME LEVEL n
TIME LEVEL n +1
t-tn+1
TIME LEVEL n +2
t=tn+2
NODE i-1
NODE i
NODE i
I
KEY:
IMPLICIT
EXPLICIT
Figure 4-2. Finite difference method temporal discretization for node i
at time level n + 1 showing explicit (dash) and implicit
(solid) relationships.
26
-------
Both sides of (4-6) contain terms which must be evaluated at the current
timestep, tn+1, and depend on the solution ^n+1, thus (4-6) is implicit in
^n+1. Another standard time differencing scheme is the Crank-Nicolson
procedure obtained by setting a=0.5 in (4-4).
The explicit scheme (a=0) is conditionally stable (the solution will not
always converge) and should only be used with small time steps. The fully
implicit scheme (a=l) is unconditionally stable although accuracy is greatly
affected by time step size. The Crank-Nicolson scheme (ct=0.5) is also
unconditionally stable and is more accurate for many flow problems [Finder and
Gray, 1977]. Theoretically, any 0 ^ a <_ 1 is first order accurate except
a=0.5, which is second order accurate for constant time steps.
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. This equation, written in matrix notation, is
.n+a
At
n+a
™n+a
i i "
n+1 _ I -
At
- (l-a)Bn+a I hn
(4-7)
where hn is the vector of known values of potential at the last time step
and hn+l is the unknown vector to be determined. The matrix A is determined
by the grid geometry and by the storage coefficients, C( 40 . The matrix B is a
function of hydraulic conductivity and of the grid geometry. The vector f is
a forcing vector that accounts for boundary conditions. It also includes
gravity flow terms. For a one-dimensional problem, the matrix B is
tridiagonal, i.e.,
B =
0
r,,
All O's
All O's
Pn-l
0
n-l
0
(4-8)
while A is diagonal (only the diagonal elements are nonzero). Thus for a = 0,
the i'th line of (4-8) expresses an implicit relation among the new potential
at node i (through q^) , the new potential at node i-1 (through p^) , and
the new potential at node i+1 (through r^). For ct=Q, the i'th equation
gives the new potential at node i explicitly.
Equation (4-7) is not a linear equation for the unknown, hn+l, since
the storage and flux matrices, An+l and Bn+l, are themselves functions of
hn+l. In order to solve (4-7) at a given time step, therefore, a variety of
27
-------
linearization and iteration schemes have been developed. The resulting linear
matrix equation can be solved either directly or iteratively, the former being
more economical in one-dimensional applications.
4.1.3 FDM Application to Vertical Unsaturated Flow—
The FDM has been applied to vertical unsaturated flow by many authors.
Freeze [1969] investigated natural ground water recharge and discharge
mechanisms. Brutsaert [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.
4.2 Finite Element Method
The Finite Element Method (FEM) has only more recently become a common
tool in ground water flow simulation. This method is conceptually more
difficult to grasp and is practically speaking more complicated than the FDM.
Nonetheless, it has several advantages including flexibility in grid design,
and improved accuracy for many nonlinear and frontal (sharp changes) problems.
Like the FDM, the FEM also involves discretization, breaking the flow
domain into subdomains. For the FEM, these spatial subdomains are called
elements and are delineated by nodes (see Figure 4-3) . Whereas the FDM
approximates the governing equation with discrete difference terms, the FEM
uses an integral approach. First, the governing equation is converted to an
integral equation, including the boundary conditions. Next, the dependent
variable (matric potential) is approximated using the values at the nodes and
interpolation functions between nodes on the elements. This approximate
potential profile is substituted into the integral equation which can be
evaluated over each element. The resulting element integrals are combined and
a system of equations is developed relating the rate of change of matric
potential at each node to the values at all nodes and to the boundary
conditions. Although the FEM can then also be applied to the integration of
the solution in time, the FDM is most commonly applied to the time derivative,
as described above. Finder and Gray [1977] describe the FEM and its
application to unsaturated flow problems.
4.2.1 FEM Spatial Approximation—
The vertical flow domain is discretized as shown in Figure 4-3. In this
domain, the matric potential is approximated as
n
y ± £ N y. (4-9)
28
-------
GRID
; - INTERPOLATION
FUNCTIONS
TOP
OF DOMAIN'
SOIL 1
SOIL
BOUNDARY
SOIL 2
BOTTOM
OF
DOMAI N
NODE 1 .
NODE 2
A
NODE
NUMNP
ELEMENT 1
ELEMENT 2
ELEMENT NUMEL
NUMEL = number of elements
NUMNP = number of node points = NUMEL +1
Figure 4-3. Finite element method discretization and linear interpolation
functions.
29
-------
where N^ is an interpolation function (see Figure 4-3) which is a function
of space, z, and ^£ is the approximate value of pressure head at node i. A
common type of interpolation function is the linear or chapeau (hat)
function. On any element, defined by the end nodes, the pressure head then
varies linearly between these two nodes (see Figure 4-3). Thus, when all the
tjji's are known, the pressure distribution is a continuous function made up
of piece-wise linear (straight) segments. Quadratic interpolation, an
alernative approach, uses three nodes per element and results in curved
segments.
In any weighted residual method, the governing equation is multiplied by
a weighting function and the results integrated over the spatial flow domain.
The weighted residual method most commonly used in ground water and soil
physical models is that of Galerkin, which uses the interpolation functions in
(4-9) as weighting functions also. This representation of (2-1) can be
written:
I
c _ _ K - dz = 0 for each k (4-10)
az az )
k ( at az
Equation (4-10) can be integrated by parts to yield:
I
'N, C -^ + K —- -^ - —~ K dz
'«, v —r + K- ' M
i k at az az az )
(4-11)
z=z
- N
k
+ K
w
0 for each k
z=o
The second part of (4-11) is non-zero only at the boundary elements. If a
flux boundary condition is used, this term, which is simply the moisture flux
at the boundary, is specified. If the fixed pressure boundary condition is
used at node k, then no solution is required at that node. At an internal
node, this term can be dropped for notational simplicity:
/ ( Ol Cl I O IN
/ 'N, C -W + K — -f^ - -r-£ K } dz = 0 (4-12)
J. , \ k at az oz 3z i
internal ( '
The value of C(tjj) and K(^) can also be approximated using the interpolation
(or other) functions:
n n
CU) = Y N.C.; K(ip) = E N-K- (4-13)
T f J T -i * ^"™ T 1
j=l J J j-1 J J
where subscript j implies the value at node j. Inserting (4-9) and (4-13)
into (4-12) we obtain:
__ / ( ai.N.^. /3N. 3N.4). \)
E J N.£.H.c. -^ + E-N.K. —L-r^+ x\ d, . 0
(e) (e) ( k J J J 3t J J J^9Z 9Z /'
/,
(4'
where the summation is over all elements.
30
-------
There are n equations, indexed by k, for n unknowns, indexed by i, in
(4-14). The time integration of (4-14) is treated in the same way as already
discussed for the FDM. The result is a nonlinear system for algebraic
equations,
'
„ « hn+ - - U-cO B hn + f na (4-15)
in which hn+l is the vector of unknown values of matrix potential. This
equation and its solution are perfectly analogous to (4-7) for the FDM. The
most significant difference is that A' , unlike A, is tridiagonal. This means
that (4-15) never gives an explicit set of equations for the new pressures,
even when a is zero.
4.2.2 FEM Applications to Vertical Unsaturated Flow —
Finite elements have been applied to vertical unsaturated flow problems
by many authors, including Neuman [1973], Johnson et al. [1982], Yeh and Ward
[1980], Hamilton et al. [1981], and Segol [1976]. Maslia and Johnston [1982]
used a two-dimensional saturated-unsaturated finite element model to
investigate ground water flow beneath a failed landfill. Trautwein, et al.
[1982] investigated leakage from a waste pond using a one-dimensional finite
element model of unsaturated vertical flow. The FEM for vertical unsaturated
flow has been verified by Segol [1976] and Milly [1982], who showed agreement
of FEM results with Philips [1958, 1969] quasi-analytical infiltration
solution.
4.3 Discretization of Flow Domain
The first step of modeling is to discretize the flow domain, in this case
the vertical coordinate and time. In general, the accuracy of numerical
techniques increase with smaller subdomains. A consistent technique converges
to the exact solution as sub domain size becomes infinitely small [Finder and
Gray, 1977]. By varying the size of these subdomains to suit a particular
simulation, accuracy can be maximized while retaining large subdomains in less
critical areas.
4.3.1 Grid Design —
The simplest grid divides the liner and site soil into equally sized
subdomains. 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
(see Figure 2-2). As infiltration proceeds, this sharp front moves into the
31
-------
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-4 shows two graded grids of
this general design. In the first, the gradient involves blocks of segments
which become larger in steps. The second example 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. The extent of
gradation and the ease of incorporating this into the model data depend on
the particular program being used. In general, the FEM is particularly well
suited for graded grids of the second type.
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.2 Time Stepping—
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. Many computer programs incorporate this feature
directly into their integration procedures. 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.
4.4. Soil properties in numerical models
During numerical solution of the governing equation or vertical
unsaturated flow, the specific moisture capacity, C(^), and unsaturated
hydraulic conductivity, K(^) must be determined. The values must be evaluated
for each value of matric potential at the node points. There are two general
methods for incorporating these continuous functions into the numerical
model. The first is by specifying a table of points from the data curves.
The computer program then interpolates between these tabular values when
needed for additional points. The second method is to directly compute soil
properties using functional relationships such as those presented above in
sections 2.3 and 2.4.
4.4.1 Tabular Interpolation—
The continuous moisture characteristic 9(^), and hydraulic conductivity
KC1^), curves can be represented by a table of values. The specific moisture
capacity, C(^) = d9/d^, is then determined from the input e(^) curve. This
32
-------
BLOCK GRADED
CONTINUOUS GRADATION
BLOCK 1
AZ =
BLOCK 2
BLOCK 3
AZ=5AZX
SUBDOMAIN 1
SUBDOMAIN NUMEL
AZ
NUMEL
NUMEL = number of subdomoins
or elements
Figure 4-4. Two vertical grids with variable subdomain sizes.
33
-------
table contains a list of values of matric potential with the corresponding
moisture contents and conductivities. When the matric potential falls between
two of these tabular points, the moisture content or conductivity is obtained
by interpolating between the respective tabular values. The computationally
simplest interpolation procedure is linear weighting. Using this procedure,
the moisture content at a given matric potential can be written:
(4-16)
where ^m and ^m+i are two values of matric potential from the input table,
with (jjm ^_ ^ilfo+i» and 9(^m) and 6(^m+i) are corresponding values of
moisture content. The specific moisture capacity for this curve is:
6 (41 .) - e dp )
m+1 m /,
,.,
m 4>m+l-4m m+1 m
Note that this parameter, which is the storage term in the governing equation,
is constant between the tabular values and is discontinuous when ^= ifm.
This is shown in Figure 4-5. Some programs actually use tabular values of
specific moisture capacity C(^m), CG|Jm+]_) and linearly interpolate for
G(IJJ). Unfortunately, this form is not the derivative of (4-16) and mass
balance errors may occur. This error will decrease as the spacing between the
tabular values decreases.
The unsaturated hydraulic conductivity can also be determined by linear
interpolation:
i
r ~ T
K(ijj) = K(IJJ ) + — K(ip ,) - K(tp ) (4-18)
m — m+1 m \
ip m+1 tp m
where K(ifjm) and KGpm+i) are the conductivities corresponding to the input
matric potentials ^m and
Log-linear interpolation is suggested by the form of functional
relationships discussed in Section 2. The moisture content between two input
points can be determined using semi-logarithmic weighting:
- log(-4>
m
m+1 m
The corresponding specific moisture capacity for this representation is:
(4-19)
C(K<) = -. - 1 m v , m/ , v [ ^ln(lO)]" (4-20)
---
which is not constant between tabular values, but is still discontinuous as
shown in Figure 4-5. Linear (4-16) and semi-logarithmic (4-19) interpolation
for moisture characterization are compared in Figure 4-5 and Figure 4-6 for
the same input table of values. Linear interpolation produces smoother
changes in C(^) .
34
-------
oi-
1
Kftl ~~
1
o
o "
^^ 1 A
rr
? (?"
i ^-
i
o f-
o
s>
'.0
LINEAR
N V — LOG
1 X X
1 X N t
1 X •* N.
\ 1 vx x
\ v NN
X
\
^ x
X
N'
X
\
X
N
.5 1.8 1.5 2.8 2.5 3.
PF
Figure 4-5. Interpolation of specific moisture capacity from table of
values 0.5 pF apart, pF = log (-Y) , Yin cm.
35
-------
LO
rr>
CO
LA
Ol
G)
•
OJ
LA
It)
.28
LINEAR
LOG
.25
.38
—r
.35
.48
.45
.58
MOISTURE CONTENT
Figure 4-6. Interpolation of moisture content from table of values 0.5
pF apart, pF = (-40 , 4» in cm.
36
-------
Log-log interpolation is suggested for the unsaturated hydraulic
conductivity. This can be written:
KOIO = KOIO - (4-21)
where
Km+l
B = —
Linear (4-18) and log-log (4-21) interpolation for unsaturated hydraulic
conductivity are compared in Figure 4-7 using the same input table of values.
The improved representation of conductivity using log-log interpolation is
obtained at a small additional computational cost.
Higher order techniques are available for interpolation. These prove
most useful in determining the specific moisture capacity from input values of
moisture content. Hamilton [1979] used cubic spline interpolation to produce
a moisture content curve, from input tabular values, that was continuous in
both first (de/dijj = C(TJJ)) and second derivatives (d^e /diJj2=dC/dijj). Thus,
the computed specific moisture capacity is smooth and continuous. This
procedure fits the input data to a polynomial with 4(NP-1) parameters where NP
is the number of tabular values. The cost effectiveness of this method
decreases as the number of input data points increases.
4.4.2 Functional Relationships—
Any of the functional relationships presented in sections 2.3 and 2.4
could be directly incorporated into a computer program. Examples of this are
Milly and Eagleson [1980] and Yeh and Ward [1980]. Two advantages of this
method are that normally few input variables are needed and it is possible to
have continuous smooth curves. The latter is a particular advantage in
developing the specific moisture capacity curve, C(ip )=dQ Aty . The
disadvantages of this method are that it can often require much more
computation than an interpolation procedure and that the functional parameters
must be determined.
The functional relationship used must closely match the soil property
data. Not only must the parameters of the function be determined, but the
functional relationship to be used must be selected. For a given function,
the parameters are adjusted to minimize some measure of the deviation between
the soil properties predicted by the model and the observed properties. For
example, a 1east-squares fit minimizes the sum of the squared deviations:
nd
(4-22)
37
-------
Figure 4-7. Interpolation of unsaturated hydraulic conductivity (K, cm/s)
from values spaced 0.5 pF apart, pF = log(-40 , ¥ in cm.
38
-------
in which K(i|>)pred. is the conductivity predicted by the model, K(^)
-------
and by (4-23)
V(i|»> = 0(4,) 2£ (4-25)
a t
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 (4-24). The spatial domain is discretized into a mesh centered grid (see
Figure 4-la) with node points on the boundary between soil elements.
Substituting the FDM represetations for the two flux terms (see (4-1) and
(4-3)) into (4-24) we have:
i+1/2 ""• x-1/2
Az.
i
(4-26)
Az.
i
where subscript "i" designates node number and
Azi+l/2 = Zi+l * zi
AZi-l/2 = Zi -
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=
(4-28)
1 /?
X K(Vl)]
A weighted temporal difference expression is obtained by substituting
(4-25) into (4-4):
in which 0 _
-------
Equation (4-29) is a FDM expression of the governing equation (4-23) at a
node. Since both Cn+l and Kn+1- depend on the solution, ^n+1, (4-29) is
nonlinear in $ . An iterative procedure is used to solve for (j;n+l. The
solution at a new iteration can be expressed as the solution from the last
iteration plus a correction computed at the new iteration
i.k+1
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 (4-29) becomes:
(4-31)
~n+a|
-,k,n+a(
(4-32)
The second term on the right hand side (RHS) of (4-29) is known from the last
time step. The first RHS term is evaluated by substituting (4-31) into (4-26):
i+l
V . V V ) -
L"i+l/2\ AZ.+1/Z
(k , k+1 ,k . , k+l\-|
•It . + Alii . - ^ • i - Aijj . \
ri i i-l 1-1 1
AZi-l/2 /_
i ,
Az.
i
Ic k.
(4-33)
Az.
or
k+1. =
/
k A*
v \
Ki+l/2\
k+1 A k+11
— A\|J .
i+1 i
Azi+l/2 >
V. (
(4-34)
in which V(ipk) £s evaluated from the last iteration, and, again, superscript
"k" or "k+1" represents iterative values at time step n+1.
The final form of the FDM governing equation is obtained by substituting
(4-32) and (4-34) into (4-29) and grouping Atjjk+1 terms on the LHS:
rk,n+a i
Ci ~AT
lk - /n
i__l
(4-35)
41
-------
All terms on the RH.S of (4-35) are known from the last time step, "n", or the
last iteration, "k" at the new time step "n+1". Equation (4-35) is solved
directly for A
-------
NEW TIME STEP n+1
COMPUTE V(tn)
SAVE <^n
CALCULATE NEW TIME STEP
At
EXTRAPOLATE TO NEW
ijin+1 USING V(ijjn)
COMPUTE FINAL
SOIL PROPERTIES
CHECK MASS BALANCE
OUTPUT SOLUTION
IF REQUESTED
END OF
SIMULATION
PERIOD?
COMPUTE SOIL PROPERTIES
FOR ipn+1
COMPUTE V(tjjk)
BUILD STIFFNESS MATRIX
LHS OF GOVERNING EQUATION
BUILD FORCING VECTOR
RHS OF GOVERNING EQUATION
SOLVE GOVERNING EQUATION
USING THOMAS ALGORITHM
HAS
SOLUTION
CONVERGED?
Figure 4-8. SOILINER solution procedure flowchart.
43
-------
At steady state no-flow, the piezometric gradient is zero:
- 0 (4-37)
Thus, the capillary pressure head gradient is equal to -1:
-r^ = - -r5- = -1 (4-38)
oz oz
The boundary conditions for this situation are
the soil column top. Between the column bottom and top, capillary
pressure head is inversely proportional to elevation:
ip = z -z z .1 z J; z (4-41)
Equation (4-41) is an analytical solution to (4-36) with boundary conditions
(4-39) and (4-40).
For the SOILINER simulations, the top node (zfc = 0) capillary pressure
head is ty = -50 cm and the bottom node (zw =-50) capillary pressure is
ty — 0 cm. Simulations with both constant node spacing and variable node
spacing produce the exact analytical solution to five significant digits.
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-42)
where a and b are constants, and n is 1, 1/2, 2, 3, or 4. For n = 4, the
analytical solution to (4-36) is given implicitly by:
/ 2 ^Jj~+ 2
In l-rr ; r I +
\ 4*2 P^ r; P2
(4-43)
tan l l-^^r II + W
44
-------
in which
r = q/a
p4 = B/r
(3 = 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-42) for the functional soil hydraulic conductivity, SOILINER was
run with boundary conditions of
^ = 0 at zw = -50 cm
and (4-44)
ty = -200 at z = 0
The flux computed by SOILINER was then used to compute Gardner's solution
(4-43) using the computer program listed in Appendix F.
Figure 4-9 shows the agreement between SOILINER simulations with a
regular grid (Test 2A) and a graded grid (Test 2C) and Gardner's analytic
solution.
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~5 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
ty = -600 cm all z, t <_ 0
ijj = 25 cm z = 0,t>0
ip = -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-10 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-11 shows the similar accuracy of
SOILINER with a graded grid. The graded grid has fewer nodes, but because of
45
-------
5>
I
O)_
G>
l
IA.
I
TEST Eft
TEST 2C
flNflLYTIC
-ffl-
.8
.5
1.0
1.5
PF
2.0
2.5
3.0
Figure 4-9. Comparison 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(o) is a regular spaced grid with
50 nodes, Test 2C(+) is a variable spaced grid with 46 nodes.
(pF = log (-40 , ^ in cm.)
46
-------
MOISTURE CONTENT
N
Figure 4-10.
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, 10^, 4 x 104, 105, and 2 x 10^ sec.
47
-------
MOISTURE COMTEMT
25 .38 .35 .40 .45
Figure 4-11.
Comparison of Philip's quasi-analytic solution (a) and
SOILINER graded grid (46 nodes) solution (solid line) for
infiltration into Yolo light clay under ponding. Curves are
at 103, 104, 4 x lO4, 105, and 2 x 105 sec.
48
-------
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.
49
-------
5. SOIL LINER DESIGN
5.1 Introduction
The usefulness of a numerical model for assessing the adequacy of soil
liner designs can be demonstrated by analysis of a hypothetical liner
performance problem. A temporary surface impoundment for storage of hazardous
wastes is to be constructed on a sandy soil with a shallow water table
5 meters below the land surface. Site soil will be excavated to sufficient
depth to install a natural soil liner. The top of the liner will correspond
to the original land surface. The impoundment depth will be 1 meter. This
impoundment is to be used to store hazardous wastes for 5 years, at which time
the liquids are to be removed and all contaminated soils excavated. With the
design life given, an accurate estimate of the required soil liner thickness
and of the extent of release at closure is needed.
5.2 Discretization of the Flow Domain
The first step in the preparation of the computer-assisted numerical
analysis and simulation is the discretization of the spatial domain, here the
soil column through which vertical flow can occur. For the problem at hand,
the soil column from the liner top to the water table can be discretized into
120 elements as follows: the first 5 cm is divided into 10 elements, the next
25 cm is divided into 25 elements, the next 50 cm into 25 elements, the next
120 cm into 30 elements, and the final 300 cm into 30 elements. Within each
block of elements, the node spacing is constant. The computer model
automatically determines the time step size each time step in such a way that
the maximum change in pressure between time steps is about 25 cm. In
addition, the maximum allowable time step size is 11.6 days.
5.3 Soil Properties from Data and Models
The site soil is a sand with a saturated hydraulic conductivity of 9.44 x
10~3 cm/s and a porosity of 0.287 (data from Haverkamp et al., 1977). The
material for the soil liner is a heavy clay with a saturated hydraulic
conductivity of 1 x 10~7 Cm/s and a porosity of 0.495. These values are
similar to a compacted bentonite clay and are more restrictive to flow than
the Yolo light clay. The unsaturated conductivity and moisture characteristic
curves for these two soils are shown in Figure 5-1. The site soil moisture is
initially in static equilibrium with the water table. Under these conditions
the pressure head at any point is the negative of the elevation above the
water table. For example, at an elevation of 300 cm above the water table,
the pressure is ty = -300 cm. The clay liner is initially at a saturation of
about 50 percent which corresponds to a moisture content of about 0.25 and a
pressure head of -500 cm.
In actuality, soil property data will be gathered by methods referenced
in Section 3 and prepared for the model as described in Sections 4.4 and
4.5.2. For this example, the computer program uses functional relationships
for K, i(J, and C. These functions are described by Haverkamp et al. [1977].
The moisture characteristic and relative hydraulic conductivity for the clay
50
-------
a)
oi.
i
n.
i
in-
V '
a '
_i
i\-.
oa.
i
s>
I
CLAY
SAND
.5
I
1.8
1.5
PF
E.0 2.5
3.0
b)
en
cu
(XJ
_
Q. ^
CLAY
SAND
00
8 5 .10 .15 .E8 .25 .30 .35 .40 .45 .50
MOISTURE CONTENT
Figure 5-1. Soil properties of hypothetical clay (solid line) and sand
(dashed line) for liner design, a) Unsaturated hydraulic
conductivity (cm/s) as a function of capillary pressure
(pF = log(-
-------
soil are those determined for the Yolo light clay. However, the saturated
hydraulic conductivity has been reduced to represent a soil suitable for use
as a liner.
5.4 Simulation
Table 5-1 shows estimates of liner thickness required for this system.
The transit time and Green and Ampt models are simplified techniques which do
not require numerical procedures (see Cogley et. al. [1982]).
5.4.1 Infiltration—
Figures 5-2 and 5-3 show the results of simulation of moisture movement
into the liner system with a liner thickness of 180 cm. This initial
thickness was determined by application of the Green and Ampt infiltration
model with a wetting from pressure head of -32 cm. Figure 5-2 shows the
capillary pressure at several locations in the liner over time. Moisture
moves into the soil at a rate which decreases with time, reaching a depth of
45 cm in about 1/2 year but reaching a depth of 90 cm after about
1 1/2 years. Likewise, at a point 135 cm below the liner top, the soil does
not approach saturation until over 3 1/2 years. The pressure at the bottom of
the liner increases in early times because the initial pressure is lower than
the pressure in the underlying sand, thus moisture actually moves up into the
liner. The pressure at the bottom of the liner first responds to infiltration
from the impoundment after about 5 years and is approaching its steady state
value at 6 years.
Figure 5-3 is a plot of the moisture content profile in the clay liner
and in the underlying sand at several times during infiltration. Before the
wetting front reaches the bottom of the liner, moisture is flowing up into the
clay from the sand due to the initial pressure gradients. The infiltration
profile at 5 years just reaches the bottom of the liner. The profile at
6 years is essentially in steady state and will remain unchanged as long as
the boundary conditions do not change. Under these conditions, there is a
steady flow of water downward from the liner toward the water table.
5.4.2 Breakthrough—
At the end of 5 years, liquid is flowing out of the bottom of the liner
at a discharge rate of 3.78 x 10~^ cm/day or 0.54 ft3/day/acre. For an
impoundment 1 acre in area, 0.54 ft^ of liquid flows out of the liner in
1 day. Whether this rate is large enough to represent breakthrough of the
liner is not clear. At this time, the moisture content in the clay at the
liner bottom is 0.275. One year later the discharge rate at the liner bottom
has increased to 1.36 x 10~2 cm/day or 19.5 ft^/day/acre. The moisture
content in the clay at the liner bottom is 0.31. At 6 years, this system is
essentially in steady state and the moisture content at the liner bottom will
riot increase significantly. As shown in Figure 5-3, the increase in moisture
content in the sand which occurs because of this infiltration between year 5
and year 6 is insignificant. The sand is so much more conductive to flow than
the clay that only a small change in pressure gradients is sufficient to move
all leakage from the liner down to the water table.
52
-------
TABLE 5-1. COMPARISON OF SIMPLLFIEP MODELS FOR 5-YEAR LINER
Pressure parameter Liner thickness
(cm) (cm)
Transit time h, = 0 74.6
d
h, = -10 77.2
d
h. = -100 97.4
d
h. = -500 155
d
Green and Ampt ^ = -10 164
ih = -32 175
rc
4> = -100 205
53
-------
1.8
2.8 3.B
TIME YEARS
4.8
5.8
£.8
Figure 5-2. Simulated pressure versus time at several points in liner
with initial moisture content of about 0.25.
54
-------
.0
MOISTURE COMTEMT
.1 .2 .3 .4 .5
63
03
63
63
(U.
l
(U-
s>
03
m.
63
in
m.
63
CD
K>
in
63
63
11.6 days
II6 days
1.6 years
6 years
Figure 5-3. Simulated moisture content profile in liner and site soil
with initial moisture content of about 0.25.
55
-------
5.5 Adjusting Liner Specifications
The effect of initial moisture content in the soil liner can be
investigated by simply changing the initial specifications on the clay soil
liner. Figure 5-4 shows the results of simulation on the same liner system
with an initial moisture content in the liner of about 0.30, which corresponds
to a pressure head of -200 cm. The fact that pressure gradients will be less
in this case might suggest that breakthrough would occur more slowly.
However, as shown in Figure 5-4, the pressure at the liner bottom begins
changing after 4 years as opposed to 5 years above. This is because the
unsaturated hydraulic conductivity is higher at this higher initial capillary
pressure, and there is less pore space which must be filled by the advancing
moisture front. This liner reaches equilibrium after 5 years. Figure 5-4
also shows that the initial pressure gradient at the liner bottom has been
reversed from that above. In this case, the pressure in the liner is higher
and the initial flow is down into the site soil, thus reducing pressure at the
liner bottom.
Use of the numerical model allows simulation of complex liner systems and
variable boundary conditions. For this third simulation, the liner consists
of a 60 cm thick layer of the site soil between two 60 cm layers of the liner
clay. In addition, this simulation incorporates a change in the impoundment
depth from 100 to 200 cm after 2 years. Figure 5-6 shows the liner moisture
content profile at three times during infiltration. Before the moisture front
reaches the upper sand layer, it has no effect on infiltration and the curve
is identical to Figure 5-3. The profile at 1.6 years, after the moisture
front has reached the sand layer, is down to about 1.5 meters below the liner
top. This compares to less than 1 meter with a homogeneous liner. Leachate
entering the sand layer quickly moves vertically down to the interface between
the sand and clay. This interface could serve as an effective leachate
collection point. The sand layer does not completely saturate until about
3 years after construction, at which time the system is essentially in steady
state. Steady state discharge from the bottom clay liner into the underlying
site soil is higher for this liner design, with a value of about 40 ft^/day/
acre, due in part to the increased impoundment depth.
5.6 Summary
The initial liner design, a homogeneous layer 180 cm thick, resulted in
release of potentially hazardous liquids to the unsaturated zone above the
water table after 5 years. If the design goal was a liner that released zero
liquid at this time, the model could be run with a thicker liner, and the flow
rate again examined at the liner bottom. In addition, numerical models can
easily simulate different designs with various soil layers and the effect of
variable boundary conditions. Soil properties may also vary from one node to
another, thus allowing simulation of a liner with variable soil properties.
The numerical model is also valuable in predicting the steady state moisture
profile and discharge rates from the liner for any configuration.
56
-------
.0
1.0
2.0 3.0 4.0
TIME YEfiRS
6.0
Figure 5-4. Simulated pressure versus time at several points in liner with
initial moisture content of about 0.30.
57
-------
MOISTURE COMTEMT
.8
.1
.2
.3
I
.4
CLAY
116 DAYS
•—• 1.6 YEfiR
5 YEflRS
S)
03
U).
SAND
CLAY
Figure 5-5. Simulated moisture content profiles in a three-layer liner and
underlying site soil.
58
-------
6. SUMMARY
Procedures for modeling flow through clay liners have been presented as
accurate and flexible tools to assist in liner design. The conceptual model
of vertical unsaturated flow has been reviewd. This non-linear partial
differential equation has metric potential,or capillary pressure as its state
variable, and requires the soils' unsaturated hydraulic conductivity and
specific moisture capacity as functions of pressure. This conceptual model
has been verified by comparison to laboratory and field tests. This one-
dimensional model is conservative for homogeneous soils.
Laboratory and field tests to determine the required soil properties are
available and have been discussed. Finite difference and finite element
numerical procedures for solving the unsaturated flow equation have been
presented. During numerical solution, soil properties are evaluated using
tabular interpolation of functional forms. Discretization of the space and
time domain for solution has been reviewed. An example computer program is
presented and tested by comparison to analytical solutins of the unsaturated
flow equation.
In order to illustrate procedures for assessing the adequacy of a single
sort liner, an example computer program (SOILINER) was applied to a
hypothetical liner system. This system consisted of a low conductivity clay
soil underlain by pervious sand. Given the soil properties of this particular
system and a design life of 5 years, a required liner thickness of two meters
was estimated. Breakthrough can be defined as the exceedence of a particular
flux rate at the liner bottom. Additional simulations demonstrated how
numerical models can be used to simulate flow through complexly layered liner
systems. Besides evaluating time to breakthrough on the basis of flux rate
out of a given liner, the numerical model also allows prediction of the flux
steady-state moisture profile and discharge rates from any soil layer of the
system under study.
59
-------
7. REFERENCES
Bear, J., Hydraulics of Groundwater, McGraw-Hill, New York, 1979.
General text with emphasis on mathematical expressions of moisture
movement. Includes a description of finite elements applied to
unsaturated flow equation.
Bouma, J, D. I. Hillel, F. D. Hole, and C. R. Amerman, Field measurement of
unsaturated hydraulic conductivity by infiltration through artificial
crusts, Soil Sci. Soc. Amer. Proc., 35, 262-264, 1971.
Brooks, R. H. and A. T. Corey, Hydraulic properties of porous media,
Hydrology Paper No. 3, Colorado State University, Fort Collins, March
1964.
Brown, K. W. and D. C. Anderson, Effects of organic solvents on the permeabil-
ity of clay soils, draft report in fulfillment of EPA Grant # R806825010,
Cincinnati, OH, 1982.
Bruch, J. C., Jr., Finite element solutions for unsteady and unsaturated flow
in porous media, California Water Resources Center, University of
California, Contribution No. 151, 1975.
Brutsaert, W. F., A functional iteration technique for solving the Richards
equation applied to two-dimensional infiltration problems, Water
Resources Research 6(7), 1583-1596, 1971.
2-D vertical, non-steady, 5 different soil types, user-oriented,
finite-difference model. New Mexico Tech., 1971.
Clapp, R. B., and G. M. Hornberger, Empirical equations for some soil
hydraulic properties, Water Resources Research, 14, 601-604, August 1978.
Added continuous gradual air entry representation to Brooks & Corey power
law model. Use Campbells relative conductivity formula. Applied model
to Holtan et al lab data including clays.
Cogley, D. R., D. J. Goode, and C. W. Young, Review of transit time equation
for estimating storage impoundment bottom liner thickness, final report
in fulfillment of EPA Contract No. 68-02-3168. GCA/Technology Division,
Bedford, MA, 1982.
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.
Finite difference on 2-D unsaturated equations and application to well
hydraulics.
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.
60
-------
Davidson, J. M., L. R. Stone, D. R. Nielsen, and M. E. Larue, Field measurement
and use of soil-water properties, Water Resources Research, 5, 1312-1321,
December 1969.
Elzeftawy, A. and K. Cartwright, Evaluating the saturated and unsaturated
hydraulic conductivity of soils, in T. F. Zimmie and C. 0. 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.
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, February 1969.
Freeze, R. A. and J. 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, 85, 228-232, 1958.
Gardner, W.R ., Field measurement of soil water diffusivity, Soil Sci. Soc.
Amer. Proc., 34, 832-833, 1970.
Giesel, W., M. Renger, and 0. Strebel, Numerical treatment of the unsaturated
water flow equation: comparison of experimental and computed results,
Water Resources Research, 9, 174-177, February 1973.
Finite difference technique on pressure and experiment on sand column,
one dimensional.
Gillham, R. W., A. Klute, and D. F. Heermann, Hydraulic properties of a porous
medium: measurement and empirical representation, Soil Science Society
of America Journal, 40, 203-207, March-April 1976.
Presents an empirical extension to King's [1965] hysteretic curve fitting
model.
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.
Hamilton, J. M., Measurement of permeability of partially saturated soils,
M. S. Thesis. Geotech Eng. Rep. GE 79-4, Univ. of Texas at Austin, 1979.
Hamilton, J. M., D. E. Daniel, and R. E. Olson, Measurement of hydraulic
conductivity of partially saturated soils, in T. F. Zimmie and C. 0.
Riggs, eds., Permeability and Groundwater Contaminant Transport, American
Society for Testing and Materials, 182-196, 1981.
61
-------
Haverkamp, R., M. Vauclin, J. Touma, 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.
Haxo, H. E., Jr., Evaluation of selected liners when exposed to hazardous
wastes, 102-111, in Residual Management by Land Disposal, Proceedings of
the Hazardous Waste Research Symposium, Tuscon, AZ, EPA-600/9-76-015,
U.S. EPA, Cincinnati, OH, 1976.
Hillel, D. I., Soil and Water - Physical Principles and Processes, Academic
Press, New York, 1971.
General text including description of unsaturated flow and moisture
characteristics. Describes finite difference techniques applied to flow
equation.
Hillel, D. I., Fundamentals of Soil Physics, Academic Press, N.Y., 1980.
Hillel, D. I., and W. R. Gardner, Measurement of unsaturated conductivity and
diffusivity by infiltration through an impeding layer, Soil Sci., 109,
149-153, 1970.
Jackson, R. D., R. J. Reginato, and C. H. M. Van Bavel, Comparison of measured
and calculated hydraulic conductivites of unsaturated soils, Water
Resources Research, 1(3), 375-380, 1965.
Validated Child, Collis-George, Marshall, and Millington and Quirk
models. Sand and loams.
Jeppson, R. W., W. J. Rawls, W. R. Hamon, and D. L. Schreiber, Use of
axisymmetric infiltration model and field data to determine hydraulic
properties of soils, Water Resources Research, 11, 127-138, February 1975.
Johnson, T. M., S. A. Rojstaczer, and K. Cartwright, Modeling of moisture
movement through layered covers designed to limit infiltration at
low-level radioactive waste disposal sites, Presented at American
Geophysical Union Spring Meeting, Philadelphia, May 31-June 4, 1982.
King, L. G., Description of soil characteristics for partially saturated flow,
Soil Science Society of America Proceedings, 29(4), 359-362, 1965.
Klute, A., Laboratory measurements of hydraulic conductivity of unsaturated
soil, 253-261, in C. A. Black (ed.), Method of Soil Analysis, Amer. Soc.
of Agronomy, Madison, WI, 1965.
Kunze, R. J. and D. R. Nielsen, Finite-difference solutions of the
infiltration equation, Soil Science, 134(2), 81-88, August 1982.
Maslia, M. L., and R. H. Johnston, Simulation of ground-water flow in the
vicinity of Hyde Park Landfill, Niagara Falls, New York, U.S.G.S.
Open-file Report No. OF82-0159, April 1982.
62
-------
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.
Analysed moisture retention data from literature and developed a
log-linear representation pF vs 0. This technique allows approximation
of entire curve from sparse data.
Milly, P. C. D., Moisture and heat transport in hysteretic, inhomogeneous
porous media: a matric head-based formulation and a numerical model,
Water Resources Research, 18(3), 489-498, June 1982.
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.
Develops 1-D vertical finite element model and verifies by comparison to
infiltration solutions. Good review of current soil moisture research.
Morel-Seytoux, H. J., Two-phase flows in porous media, 119-202, in V.T.
Chon, 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.
Integral expression using soil moisture data, comparison of technique to
others (Childs, Collis-George, Averjanor, Gardner, Millington and Quirk)
for 45 soils.
Mualem, Y., Hydraulic conductivity of unsaturated porous media: generalized
macroscopic approach, 14(2), 325-334, April 1978.
Functional power law between Kr and Se using work integral from soil
moisture curve, comparison with 50 soils.
Neuman, S. P., Saturated-unsaturated seepage by finite elements, Journal of the
Hydraulics Division, ASCE, 99(HY12), 2233-2250, December 1973.
UNSAT II. Computes hydraulic heads, pressure heads, water content,
boundary fluxes and internal sinks and sources in a
saturated/unsaturated, nonuniform, anisotropic, porous medium under
nonsteady state conditions.
Nielsen, D. R., D. Kirkham, and W. R. van Wijk, Diffusion equation calculations
of field soil water infiltration profiles, Soil Science Society of
America Proceedings, 25, 165-168, 1961.
Olson, R. E. and D. E. Daniel, Measurement of the hydraulic conductivity of
fine-grained soils, in T. F. Zimmie and C. 0. Riggs, eds., Permeability
and Groundwater Contaminant Transport, American Society for Testing and
Materials, 18-64, 1981.
63
-------
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, New York, 1969.
Finder, G. F. and W. G. Gray, Finite Element Simulation in Surface and Sub-
surface Hydrology, Academic Press, New York, 1977.
Advanced text including description of finite element application to
unsaturated flow and examples- Good comparison between finite difference
and finite element techniques.
Ragab, R., J. Feyen, and D. Hillel, Comparison of experimental and simulated
infiltration profiles in sand, Soil Science, 133(1), 61-64, January 1982.
Compared experimental and simulated infiltration when conductivity data
obtained from water desorption method and mercury intrusion method. Lab
tests. Sand.
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.
Roberts, D. W., Soil properties, classification, and hydraulic conductivity
testing, draft report in fulfillment of EPA Contract No. 68-03-3058,
Cincinnati, OH, 1982.
Rogowski, A. S., Watershed physics: model of the soil moisture
characteristic, Water Resources Research, 7(6), 1575-1582, December 1971.
Rose, C. W., and A. Krishnan, A method of determining hydraulic conductivity
characteristics for non-swelling soils in situ and of calculating
evaporation from bare soil, Soil Sci., 103, 369-373, 1967.
Segol, G., A three-dimensional galerkin-finite element model for the analysis
of contaminant transport in saturated-unsaturated porous media, Finite
Elements in Water Resources, Proceedings of 1st International Conference,
Pentech Press, 1976.
Siegel, J. and D. B. Stephens, Numerical simulation of seepage beneath
lined ponds, 219-232, in Symposium on Uranium Mill Tailings Management,
Proceedings, Colorado State University, Fort Collins, Co., 1980.
Swartzendruber, D., The flow of water in unsaturated flows, 215-292,
in R. J. M. DeWiest, ed., Flow Through Porous Media, Academic Press, New
York, 1969.
Trautwein, S. J., D. E. Daniel, and M. W. Cooper. A case history study of
water flow through unsaturated soils, Presented at American Geophysical
Union Spring Meeting, Philadelphia, May 31-June 4, 1982.
64
-------
Watson, K. K., An instantaneous profile method for determining the hydraulic
conductivity of unsaturated porous materials, Water Resources Res., 2,
709-715, 1966.
Weeks, L. V., and S. J. Richards, Soil water properties computed from transient
flow data, Soil Sci. Soc. Amer. Proc., 31, 721-725, 1967.
Yen, G. T. and D. S. Ward, FEMWATER: a finite-element model of water flow
through saturated-unsaturated porous media, Oak Ridge National
Laboratories, Report No. ORNL-5567, October 1980.
User's manual and modifications to 2-D vertical flow model of Reeves and
Duguid (1975). Based on pressure head.
Youngs, E. G. , Moisture profiles during vertical infiltration, Soil Sci.,
84(4), 283-290, 1957.
65
-------
APPENDIX A
REVIEW OF THE TRANSIT TIME EQUATION FOR ESTIMATING
STORAGE IMPOUNDMENT BOTTOM LINER THICKNESS
This report was prepared for the
Office of Solid Waste
under contract no. 68-02-3168
A-l
-------
DISCLAIMER
This report was prepared by David R. Cogley, Daniel J. Goode, and
Charles W. Young of GCA Corporation, Technology Division, Bedford,
Massachusetts under contract no. 68-02-3168. 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.
ii
-------
CONTENTS
Figures iv
1. Introduction 1
2. Estimation of Liner Thickness 2
Derivation, assumptions, and criteria for the transit
time equation 2
Limitations of the transit time equation 4
Possible improvements or modifications to the transit
time equation 5
Green-Ampt wetting front model 8
Transient linearized approximate solution 10
Model comparisons 12
3. Unsteady Flow in Clay Soils 15
Suction properties of soils 15
Numerical solution of unsaturated flow 25
4. Summary 28
Conclusions 28
Recommendations 29
Liner specifications 29
5. References 31
ILL
-------
FIGURES
Number Page
1 Definition sketch for proposed EPA equation 7
2 Green-Ampt Infiltration Model 9
3 Graphical solution to linearized infiltration .... 11
4 Hydraulic properties of Yolo light clay . 12
5 Comparison of liner thickness equations 13
6 Method of measuring soil suction 16
7 Relationships between suction and moisture content for a silty
sand 16
8 Suction/moisture content and shrinkage relationships for a
heavy clay soil 18
9 Response of a clay/sand interface to wetting 20
10 Wetting behavior at clay/sand interface 21
11 Schematic cross section through a permeater 22
12 Comparison of measured suctions with values calculated using
measured properties and a finite-element computer program . . 22
13 Relationship between water content and suction for Goose Lake
clay 23
14 Calculated coefficients of hydraulic conductivity for Goose
Lake clay 24
-------
SECTION 1
INTRODUCTION
The Environmental Protection Agency is incorporating an alternative
design concept into proposed regulations for hazardous waste disposal
(Part 264, Subpart K) to allow construction of liquid storage impoundments
with a single natural soil bottom liner of sufficient thickness to prevent
liquid breakthrough during the operating life of the impoundment. At closure,
liquid waste is to be removed along with the depth of bottom liner
contaminated as a result of liquid waste seepage. The goal of incorporating
this concept into the regulations is to provide flexibility to the prospective
owner/operator in selecting a means of storing hazardous wastes prior to
ultimate disposal. The Agency is considering the use of a transit time
equation to provide a simple method of estimating necessary bottom liner
thickness as a function of the design impoundment life.
This report evaluates the correctness of the transit time equation being
considered by EPA for estimating liner thickness and is intended for insertion
into the Administrative Record of Rulemaking. The review identifies the
derivation of the equation, key assumptions and other criteria, applicability,
reliability, inherent limitations, and possible modifications or
improvements. Other models and experimental and engineering methods that can
be used to estimate liner thickness are presented and compared. The complex
and highly variable adsorption, molecular diffusion, and reactive properties
of individual components and their effects on liner thickness have not been
included in this review. Rather, the review concentrates on flow of a single
fluid through the liner. Conclusions and recommendations are presented to
illustrate regulatory needs.
The transit time equation under consideration is derived from Uarcy's
equation for one dimensional, steady state, saturated flow. Other basic
assumptions include the use of total (versus effective) porosity, and constant
hydraulic conductivity independent of moisture content. Thus, the intent is
to estimate liner thickness using documented or measured values for hydraulic
conductivity and soil porosity.
-------
SECTION 2
ESTIMATION OF LINER THICKNESS
The derivation, applicability, and limitations of the transit time
equation are highlighted in this section. Modifications to the expression are
discussed including procedures to incorporate effective porosity and negative
pore pressure at the liner bottom. A review of modeling efforts adopted by
Pope-Reid Associates based on the Green-Ampt equation for movement of the
wetting front through clays is then presented. An alternative mathematical
model which attempts to address unsteady state conditions (using a linearized
approximate solution) is presented. The merits and limitations of each of
these alternatives are explained in this section. Based on this review and
knowledge of the mechanics of liner wetting and saturation, we present in
Section 3 a proposed approach for measuring liner breakthrough times in the
laboratory and field and estimating liner breakthrough times by scaling up
observed behavior.
DERIVATION, ASSUMPTIONS, AND CRITERIA FOR THE TRANSIT TIME EQUATION
The transit time equation under consideration by EPA was suggested by CMA
in the development of guidance for the Part 267 permitting standards (A. Day
correspondence, June 1982). The equation takes the form:
(1)
d = 0.5 — , i-
n ( Vn / \ n
where:
d = necessary thickness of soil (feet)
n = total porosity
K = hydraulic conductivity (ft/yr) (which is a function of the soil medium
and fluid flowing through it)
h = maximum fluid head on the liner (feet)
t = facility life from startup through closure (years)
and is derived based on the illustration and reasoning presented below:
-------
WASTE
k,n
NATURAL
SOIL
LINER
WATER TABLE
The transit time equation (1) is derived from Darcy's Law, which states:
K .
v = i
n
(2)
where:
v is the velocity and
i is the hydraulic gradient
If t is the desired containment time, then the liner thickness, d, is given by
Kit
d = v • t =
n
(3)
The hydraulic gradient, i, is given by
h + d
i =
d'
, yielding
Rearrangement to set the expression equal to zero gives:
,2 Kdt Kht _ .
(4)
(5)
(6)
n
n
which can be solved for d using the positive solution of the binomial
equation, or
d = 0.5
n
(7)
-------
Key assumptions and criteria implicit in using the transit time equation
are:
• saturated flow consistent with Darcy's Law
• one dimensional vertical liquid flow
• steady state conditions
• capillary tension (suction) = 0
• pore fluid pressure at bottom of liner is equal to atmospheric
pressure
• mass transport by advection only (i.e., no dispersion or diffusion)
LIMITATIONS OF THE TRANSIT TIME EQUATION
The scenario under consideration requires provision of a single soil
liner thick enough to contain all leachate within the basin and liner
throughout the life of the storage impoundment. Therefore, an equation is
required that can estimate the time to breakthrough of leachate at the liner
bottom.
To accurately and reliably estimate the time to breakthrough for a given
liner thickness, an equation is required that accounts for variable liner
moisture content during wetting and associated variability in leachate
conductivity and suction forces induced by capillary tension at any moisture
content below complete saturation. Moore (SW-869, Sept. 1980) points out that:
"(1) during the early stages of wetting of a compacted clay liner,
capillary attracton forces will predominate over gravitational
forces; and
(2) as the clay liner becomes wetter the capillary forces decrease in
importance; and, when the liner is saturated, these forces become
negligible in comparison to gravitational forces."
The transit time equation under consideration is applicable (although
limited) in describing leachate movement after wet up but is not appropriate
to describe the first stage of liner wetting under capillary forces because:
• steady state conditions are assumed
• capillary tension is ignored
• conductivity is assumed constant independent of liner moisture
content
The use of the saturated conductivity value, which is greater than
unsaturated conductivity, results in higher values of liner thickness, and
-------
partially offsets the neglect of capillary forces. As shown at the end of
this section, this offset is not sufficient to render accurate predictions of
required thickness (see Figure 5).
Although more applicable to the phase of leachate movement after initial
wet up, additional limitations of the transit time equation must be recognized:
• total, rather than effective porosity is used
• liner homogeneity is assumed, thus ignoring localized failures such
as liner cracking
• potential changes in liner conductivity or permeability as a
function of long term waste liquid/liner interactions are not
accounted for
The net effect of these limitations is certainly nonconservative, in that
the resulting liner will be of insufficient thickness to contain leachate
(avoid breakthrough) during the operating life of the storage impoundment.
POSSIBLE IMPROVEMENTS OR MODIFICATIONS TO THE TRANSIT TIME EQUATION
Within the general framework of a steady-state advective transport
analysis, modifications to the proposed equation and its application could
increase its usefulness. These modifications are not an attempt to change the
basic assumptions which are:
• steady state Darcy flow
• fully saturated flow
• advective transport only
• homogeneous liner
The proposed liner thickness equation is modified to include the effects
of:
• effective porosity instead of total porosity
• negative fluid pressure at liner bottom
The steady state saturated vertical flow equation in the liner can be
written (see Bear 1979):
2--Hc|*K-0 (8)
9z dz
where K [L/Tj is the vertical saturated hydraulic conductivity, <}>=£+ z [L]
Y
is the piezometric head, p [F/L2] is the fluid pore pressure, Y[F/L3] is
-------
the fluid's specific weight, and z [LJ is the vertical Cartesian coordinate,
positive upward. Hydraulic conductivity is a function of both the porous
medium matrix and the pore fluid. Assuming K is constant over z, equation
(8), and Darcy's law give:
K 1& = constant = q (9)
dZ
where q [L/Tj is specific discharge, or volume flux per unit area of the
porous media. The fluid velocity is obtained by dividing the specific
discharge by the ratio of active fluid flow area to total area. This ratio is
equivalent to the effective porosity, ne, which is less than total porosity,
n, because of closed pores, dead-end pores and entrapped air pores which do
not contribute to the flow area. Effective porosity is always less than or
equal to total porosity. Thus, fluid velocity V [L/T], is
v-i- £- -I* do)
n n dz
e e
The fluid pressure boundary conditions (see Figure 1) on the liner can be
used to evaluate (10). The piezometric head at the top of the liner is equal
to the elevation of the impounded water free surface hs[L]. At the liner
bottom, the boundary condition is more ambiguous. If the fluid pressure is
assumed to be in static equilibrium with the atmosphere (pj, = 0), then the
piezometric head is equal to the elevation of the liner bottom:
where subscript b, ( )jj> implies liner bottom value.
An alternative boundary condition can be formulated based on capillary
forces. Capillary rise is the phenomenon of a fully saturated zone above the
p=0 horizon, where the pore pressure is less than zero (see Bear 1979).
Likewise, at the liner bottom, the liner may be fully saturated and have
negative pressures due to drainage to the underlying soil or vapor transport
from the interface. In terms of pressure head (p/Y; this force is called, for
example, critical capillary head, hcc (Bear 1979) or displacement pressure
head, h = h. + z, (12)
Dab
where h, h -h,-z,
ii _ 9t - Yb _ s d b (13)
3 z z-z z-z
-------
DATUM
1
h
i _
1 =r
IMPOUNDED
LIQUID
1 1
i
Z
/////x
/ LINER /
1 MATERIAL '
/ }//
/i / / /i/
1 0% 4 0 0. , 0 0 | 0 0
T o °* * 00
%• -4
°0 NATURAL b&
Zb FILL °
t
h
1
USING(4)
Zb + bd Zb
PIEZOMETRIC HEAD
Figure 1. Definition sketch for proposed EPA Equation.
-------
where zt [L] is the elevation of the top of the liner. For no capillary
forces (p = 0), h
Rearranging and solving the binominal equation for d gives:
d = 1/2
ne
+ 4
(^)
(17)
This expression is identical to the unmodified equation (7) if the fluid
pressure at the liner bottom is zero (h(j = 0) and effective porosity (ne)
is replaced by total porosity, n.
These two modifications,
• effective porosity instead of total porosity
• negative fluid pressure at liner bottom
both increase the design thickness, d, for a given design life.
GREEN-AMPT WETTING FRONT MODEL
Green and Ampt (1911) derived a simple model of infiltration which has
been proposed as a design model for liner reliability (Pope Reid Associates
Correspondense, 1982). As shown in Figure 2, the soil moisture profile is
conceptualized as a square wave moving down the soil column. Above the
wetting front, the soil is fully saturated, while below the wetting front the
moisture content is equal to its initial value, 6^. Assuming the pressure
head at the front is a suction, ^c, [L] due to the partial initial
saturation, Darcy flux q [L^/T] in the saturated zone is: (see McWhorter
and Nelson 1979).
-------
LINER TOP
LINER BOTTOM
GREEN-AMPT
UJ
_)
u
•GREEN-AMPT
-*,
MOISTURE CONTENT
SUCTION -I//
Figure 2. Green-Ampt Infiltration Model.
-------
(h+L -^ )
q = K LP C (18)
P
where Lp [L] is the depth of penetration of the wetting front from the liner
top. Conservation of volume of the pore fluid requires:
dL
q = (n-6.) —£ (19)
Combining (18) and (19) and integrating:
n-e. r
* " T-M LP - (h-V ln
(20)
The design liner thickness is equal to d = L.,, when t equals the design life.
/h+L -<|> \
P_c
\ h-) [-] is moisture content, and
K (fy) [L/T] is effective hydraulic conductivity. Note that z is positive
downwards from the liner top (z = 0). Equation (21) is highly nonlinear due
to the dependence of 9 and K on fy. These nonlinearities can be removed, and a
solution obtained, by rearranging:
3* _ d£ JL ir J*i . d£ dK ii=n
3t d6 3t 3z d6 di|> dz
and substituting D* = ^| K and K* = ^| 37 - II resulting in:
do do dip do
10
-------
at
+ K*
- 0
(23)
Applying boundary conditions (Jj = h at the top (z = 0), where h is the depth of
ponded liquid, and initial condition 4> = ^i where ijjj_ is the initial fluid
pore pressure head, the solution is: (Bear 1979, page 268)
(24)
where erfc is the complimentary error function (c.f. Crank 1975).
Equation (24) is presented graphically in Figure 3 with d = z. Unlike the
proposed transit time equation or the Green-Ampt model, the transient
linearized solution results in a continuous profile of soil moisture which is
physically more accurate (see "Actual" moisture content curve in Figure 2).
Thus there is no sudden increase in moisture content or reduction in suction
pressure at the liner bottom. Rather, moisture content gradually increases
from the initial value to saturation as the wetting front advances. This is
consistent with the experimental observation that moisture flows out the
bottom of the column before it is completely saturated (see Mclntyre, et al.,
1979). Since there is no explicit breakthrough, physically or mathematically,
we must define this concept in terms of relative changes in moisture content
or pressure. A common choice is 50 percent (see Figure 2) but obviously
significant leachate has already reached the liner bottom at this time. A
standard value of ('l'-^^) /(h-4^), such as 0.1, could be chosen to represent
break-through. Liner thickness is then found from Figure 3 iteratively by
assuming a thickness, d, computing D*/K*d, and finding K*t/d on the bottom
axis. Time t is then computed, and a new d assumed depending on whether t is
less than or greater than the design life.
10 50
Figure 3. Graphical solution to linearized infiltration
(after Ogata and Banks 1961).
11
-------
The major fault with this method is determining D* and K* since the terms
which define them, especially D*, vary over several orders of magnitude during
saturation. Moore (1982) has outlined a method for determining K** and D**
for the moisture content form of the linearized governing equation. (Note K**
and D** are not numerically equivalent to K* and D* due to different governing
equations.) Saturated conductivity is used for K**, and D** is determined by
curve fitting results of laboratory column tests. It has not been verified
that results of short—time small-scale lab tests can be adequately scaled to
represent field behavior, and the fitted D** is a function of time and length
(Daniel, 1982). This method does approximate the dynamic characteristics of
the infiltration event, specifically, capillary forces, and decreased wetting
front velocity with time.
MODEL COMPARISONS
The liner thickness design equations discussed above can be compared to a
numerical finite element solution of the nonlinear unsaturated flow equation
(Milly, 1982) and a quasi-analytical solution (Philip, 1958). Finite element
solution of the nonlinear unsaturated flow equation is described in
Section 3. It is considered that these two solutions are the most accurate,
> the "correct" result.
The liner material simulated is Yolo light clay. This material is not a
particularly effective liner, but has been well studied and is presented for
comparison. The soil properties are: n = 0.495, 6 j_ = 0.237, and K = 1.23 x
10 ^ cm/s. The liquid depth in the impoundment is h = 25 cm. The relative
hydraulic conductivity (K1( )/K) and moisture retention functions are plotted
in Figure 4. Breakthrough is defined by the ratio of suction pressure
reduction at the liner bottom to the difference between initial pressure and
the top boundary pressure. Results of the numerical and quasi-analytical
solutions are plotted for 10 percent and 50 percent reduction in suction at
the bottom of the liner (Figure 5) to illustrate the conservative effect of
choosing values lower than 50 percent.
in
H
o
£ •
d
(K
~ «i -
O
1.3
l.'4
is
PF MOISTJRE CONTENT
Figure 4. Hydraulic properties of Yolo light clay
(pF = log (-), in cm).
12
-------
numerical and quasi-analytical
10% after Milly 1982
transit time solutions (see text)
i .
1\ /,
Unmodified transit time equation j
Q Modified transit time equation '
O Green-Ampt (1911) infiltration model
A Linearized unsaturated flow
~ § -m y- ''—
10 15 20
LINER THICKNESS d , cm
Figure 5. Comparison of Liner Thickness Equations
13
-------
The unmodified EPA transit time equation (1) is used to calculate liner
thickness. Design depth is plotted (O) versus design life on a semi-log
scale. The modified transit time equation (17) is applied to the liner (d).
The additional parameters needed for the modified form are taken as: ne = n
= 0.495, for the top line hd = -10cm, and for the bottom line hd = -100 cm.
Results of the Green-Ampt model (o) are plotted for two values of suction
head, the top line is tyc = -10 cm and the bottom line is ^c = -100 cm.
The coefficients for the linearized unsaturated flow model are determined
by fitting to the numerical results, D* = 10~2 and K* = 3 x 105. This
model's results are plotted for one set of parameters (A).
As shown in Figure 5, solutions differ over a wide range. The EPA
transit time equation underestimates liner thickness for given design life.
The modified solution is also sensitive to h^, a parameter which is an
artifact of the approximations involved, and cannot be determined from field
data. This error would be even more pronounced for a heavy clay, with lower
conductivity. The Green-Ampt and linearized equation models are both very
sensitive to parameters which, again, are artifacts of their respective
approximations and cannot be accurately estimated from soil property data.
14
-------
SECTION 3
UNSTEADY FLOW IN CLAY SOILS
Neither steady state models nor linearized unsteady flow models are
capable of accurately predicting wetting phenomena for compacted clays. A
brief recounting of soil suction for a variety of soils is presented here
prior to a description of models which can properly account for unsteady flow
phenomena.
SUCTION PROPERTIES OF SOILS
The relationship of soil suction and moisture content is illustrated for
several soils. An apparatus for measuring soil suction is used to illustrate
the concept of soil suction. This is followed by an explanation of units
commonly used to express soil suction. Properties of a silty sand are used to
illustrate soil suction values for a simple case. Properties of heavy clay
are used to illustrate the effects of wetting, drying, and mixing on soil
suction. Next, the behavior of a clay-sand interface is noted. Finally,
requirements of a numerical model and its applicability to data from
laboratory and field tests are discussed. It is concluded that flow of
moisture out of the clay occurs not when the moisture content at the bottom of
a liner begins to increase but rather when the clay is almost saturated.
Regulations should address this saturation phenomena and not the arrival of
the leading edge of the wetting front.
There are many techniques for measuring soil suction. One method is
illustrated in Figure 6. A sample of soil is placed on a porous plate which
is in contact with a water filled reservoir. The water is under a
controllable vacuum (suction). At zero vacuum, the water meniscus moves
toward the soil. As vacuum is applied, movement of the meniscus slows. When
the vacuum equals the soil suction, the meniscus become stationary. This
standard technique is described in Part 19 of ASTM Standards (D2325, D3152).
Several scales or units of suction are in common usage. Soil suction may
be expressed in centimeters of water. A pF scale has been defined such that
pF equals the logarithm to the base 10 of the soil suction in centimeters (of
water). Soil suction values are also frequently expressed in atmospheres.
The following values of soil suction are equivalent: 1 atmosphere, 14.7
pounds per square inch, 1,033 centimeters of water, pF 3.0
The relationship of soil suction to moisture content for a silty sand is
illustrated in Figure 7. The drying curve is shown with open circles and the
15
-------
Sample chamber
Sample
Porous plate
Water-filled reiervoir
Vacuum lint
Water menucui''
Manometer,
Figure 6. Method of measuring soil suction
(Croney and Coleman).
5 IO 15 20 25
MOISTURE CONTENT - percent
Figure 7. Relationships between suction and moisture content for
a silty sand (Croney and Coleman).
16
-------
wetting curve with filled circles. The observed hysteresis is typical of all
soils. The important point is that for a variation of 1 pF unit in soil
suction anywhere from 2 percent to 20 percent changes in moisture content (100
times weight of water divided by weight of dry soil) are observed. The
greatest changes in moisture content occur for pF values in the vicinity of pF
1, i.e., 10 cm of water suction. For this particular sample of silty sand,
very little water will enter the sand from adjacent soils if their suction
values exceed pF 3.
Suction values of a heavy clay soil (plastic limit 27 percent, liquid
limit 77 percent) subjected to a variety of experimental conditions were
investigated by Croney and Coleman and their data are presented in Figure 8.
Their explanation is excerpted below.
"[Figure 8] shows the suction/moisture content relationships for a heavy
clay soil (London clay). Curve A represents the "undisturbed" material
drying from a suction close to zero. Inset is shown the shrinkage curve
for the soil drying from the field-moisture condition... If the drying
process is continued to oven-dryness and the suction is then decreased,
the wetting condition of the soil over the range pF 1 - pF 7 is
represented by the curve B [Figure 8J. On wetting to pF 1 the moisture
content of the soil is appreciably lower than its value prior to
oven-drying, showing that some irreversible structural change has
occurred as a result of oven drying. On drying for a second time to pF 7
the suction relationship follows curve C, which forms a closed loop with
curve B. Any subsequent wetting and drying cycles over the range pF 1 -
pF 7 give suction curves corresponding to this loop, which appears
therefore to be unique for the soil. When the structure of the clay is
partially destroyed by thoroughly slurrying the soil at a high moisture
content, and the slurried clay is subjected to an increasing suction,
pF 0 - pF 7, curve D is obtained... If the slurried soil is dried to a
suction less than pF 4.8 and the suction is then reduced progressively
and subsequently again increased, closed hysteresis loops of the type E
and F are obtained...
"During the last few years several workers have attempted to ascribe pF
values to the Atterberg limits... The heavy clay soil referred to
[Figure 8] was thoroughly slurried with water to a moisture content well
above the liquid limit. The number of blows required to close the
groove, the suction of the soil immediately after the test, and the
moisture content were all determined... Before each test the soil was
thoroughly mixed in accordance with the standard procedure. The
suction/moisture content points for various numbers of blows are shown
[Figure 8]. The suction of the plastic limit was also determined... The
liquid and plastic limit points were found to lie on the continuous
dotted line, G. It was subsequently found that if the soil was disturbed
without change of moisture content irrespective of its initial suction,
it assumed a suction given by curve G at that moisture content... The
broad conclusion is that if the suction/moisture content relationship of
the soil is represented by a point above the line G, disturbance causes a
17
-------
00
ft 12 16
MOISTURE CONTENT
Number of blowt M
35 I limit t«it
4O SO 6O 7O
MOISTURE CONTENT, percent
Figure 8, Suction/moisture content and shrinkage relationships for a heavy clay soil,
(Croney and Coleman)
-------
decrease in suction to the value given by the line at the moisture
content of the sample. On the other hand, if the initial suction is
represented by a point below line G, disturbance will be accompanied by
an increase of suction. Tests on other heavy clays have shown that for
each soil there is a unique suction/moisture content relationship."
The response of a clay/sand interface to wetting is illustrated in
Figure 9 which is derived from curve B of Figure 8 and the wetting limb of
Figure 7. Curve B of Figure 8 approximates the wetting curve for a heavy clay
placed near optimum moisture content which is expected to be a few percent dry
of the plastic limit. As discussed above, the sand will not absorb
appreciable moisture until the suction aproaches pF 1. Thus, as the wetting
front approaches the clay/sand interface, the suction will fall from high
values toward pF 0. As the clay reaches pF 5, its moisture content is 10
percent and that of the sand is less than 2 percent. At pF 3, the moisture
content is 23 percent for the clay and 2 percent for the sand. Between pF 1.5
and pF 0.5 the sand water content increases by 20 percent to 24 percent.
Somewhere in this range gravity flow of water downward through the sand is
expected. Thus, in the present case, flow of water through the sand is not
expected until water content in the clay has risen from its initial value to
28 percent. See Figure 10 for an alternative display of this phenomenon.
For a waste impoundment with a clay soil liner underlain by unsaturated
sand, it is important to determine, for the sand, the moisture content at
which free drainage occurs and the corresponding pF value. Then, if one can
model the wetting behavior of the clay, it should be possible to predict the
breakthrough time at the sand/clay interface if the initial moisture content
of the clay is known. A suitable model is available and is described in the
next subsection. It is instructive to consider the theoretical basis for the
model and then to determine how practical use may be made of the model using a
limited number of laboratory and field tests.
Hamilton, Daniel and Olson (1981) describe the successful application of
a Galerkin finite element model to the quantitative prediction of soil suction
versus depth and time as moisture moves into a clay sample previously prepared
at a low moisture content. Data were reported for compacted samples of a fire
clay known commercially as Goose Lake clay (plastic limit 19 percent, liquid
limit 26 percent) composed of 19 percent sand, 62 percent silt, and 19 percent
clay. The standard maximum dry density and optimum water content were
122 lb/ft^ and 10.5 percent, respectively. The test cell is shown in
Figure 11. Moisture moving into the clay produced changes in soil suction
which were measured by the sensors (thermocouple psychrometers).
In Figure 12 measured values of suction are contrasted with values
computed by the Galerkin finite element model. Agreement is excellent for all
but the first sensor probe. This same computer model could be used to predict
clay soil liner lifetimes based on pF values corresponding to breakthrough at
a clay/sand interface. Input data required for the model are shown in
Figures 13 and 14.
19
-------
10 20
MOISTURE CONTENT, p.rcint
Figure 9. Response of a clay/sand interface to. wetting
(after Croney and Coleman).
20
-------
GRAVITY FLOW OF WATER
IN SAND EXPECTED
SANO MOISTURE CONTENT, ptrctnt
Figure 10. Wetting behavior at clay/sand interface.
21
-------
THERMOCOUPLE
PSVCHROMETERS
O-RING
WATER
SOURCE
HYPODERMIC
SYRINGE
NEEDLE
.VENT
Figure 11. Schematic cross section through a permeameter
(Hamilton, et al.).
100 200 300
TIME, MRS
400
500
Figure 12. Comparison of measured suctions with values calculated using
measured properties and a finite-element computer program
(Hamilton, et al.).
22
-------
12
O)
o 10
0>
Q.
z
UJ
z
o
tr
UJ
I
8
T 1 1 1—» 1 i I
o
—I 1 1 fill
o WATER CONTENT TEST
• PERMEABILITY TEST
— USED IN COMPUTER
j ' ' 1—I—I I 1 I
1—I I
100
SUCTION, ATM
Figure 13. Relationship between water content and suction for Goose Lake clay,
(Hamilton et al).
-------
o
UJ
CO
o
>
H-
>
h-
u
Q
Z
O 10"'
10
,-12
• PROBE I
A ?
o 3
o 4
* 5
— USED IN
COMPUTER
20
40 60
SUCTION, ATM
80
100
Figure 14, Calculated coefficients of hydraulic conductivity for Goose Lake
clay,(Hamilton et al.)
24
-------
For a waste impoundment, a clay liner would be placed just wet of the
optimum moisture content which for Goose Lake clay is 10.5 percent. The
expected suction value is under five atmospheres (pF 3.7). Thus, far less
experimental data would be required than Figures 13 and 14 might otherwise
indicate. For the purposes of predicting liner lifetimes, it is not yet clear
whether one should gather basic data such as that shown in Figures 13 and 14,
or whether one should employ a curve fitting procedure to data of the type
shown in Figure 12.
NUMERICAL SOLUTION OF UNSATURATED FLOW
Due to the nonlinearities of the unsaturated flow equation, exact
analytical solutions have not been obtained. Thus, the approximate treatments
have been discussed above. Numerical techniques are available to solve the
unsaturated flow equations. These techniques, namely finite difference and
finite elements, retain the nonlinearity of the liner column matrix with
piece-wise continuous approximate representations and are thus most faithful
to the original governing equation. Available programs range from basic flow
(Bruch, 1975, Johnson, et al., 1982) to coupled heat and moisture transport
with vapor transport and hysteresis (Milly, 1982). Using these programs, it
is possible to evaluate specific liner designs, with parameters from field and
laboratory tests, as well as develop generic "empirical" relations between
field and lab parameters and flow characteristics. Not the least important of
these characteristics is breakthrough time for a wetting front.
The finite element method is an advanced numerical technique related to
finite difference techniques. We start with the governing one-dimensional
flow equation on our domain, the liner:
3* d9 3 3* dK 3* _ n m
9t d^ 9z "H d^ "al ~ U UJ
where ^ [Lj is the suction head; Q[l?/l?\ is moisture content; K = K (, we replace it by an
approximation
n
Ni
where N£ is an interpolation function which is a function of space, z, and
$^ are values of approximate pressure head at node i. A common type of
interpolation function is the linear or chapeau (hat) function. On any
element, defined by the nodes, the pressure head at any point is a linear
interpolation of the pressure head at these two nodes. Thus, when all the
ipi's are known, the pressure distribution is a piece-wise continuous
function made up of linear (straight) segments. Quadratic interpolation uses
three nodes per element and results in curved segments.
25
-------
The terms dQ/d4> and K (^), which depend on
-------
In addition to basic soil data (porosity, density, etc.) the numerical
model requires definition of the hydraulic conductivity and moisture retention
curves. Moisture retention is easily measured by equilibrating a soil sample
with a controlled suction and measuring water content by weight. The
hydraulic conductivity relationship is more ellusive, but can be effectively
determined (see Hamilton, et al., 1981; Daniel, 1982).
There are many methods available for generating hydraulic conductivity
curves (Miller and Bresler, 1977; Brutsaert, 1979; Clapp and Hornberger, 1978;
Elzeftawy and Cartwright, 1981; Gruber, 1982; Brooks and Corey, 1966). The
method of Elzeftawy and Cartwright (1981) involves computing the hydraulic
conductivity curve incrementally using the measured moisture retention at a
given suction. These methods have been successful with granular soils but
their applicability to clays is tenuous.
Although it is beyond the scope of this work to recommend laboratory
procedures, the evaluation of a liner can be outlined. Laboratory tests on
the liner material generate preliminary soil data, including hydraulic
properties. As the liner is installed, undisturbed samples from each lift are
analyzed to verify anticipated properties (moisture content, Proctor density,
etc.). Finally, the entire liner column is monitored and simulated using
in-place properties as determined from undisturbed samples. The model can
also be used to estimate sensitivity of the thickness to errors in parameters.
27
-------
SECTION 4
SUMMARY
CONCLUSIONS
Based on GCA's review of mathematical expressions or models applicable to
estimate the thickness of storage impoundment bottom liners, we conclude that
the transit time equation is inappropriate for predicting the time to
breakthrough. Under unsaturated conditions, capillary tension controls the
movement of the wetting front whereas the influence of liquid head above the
liner and gravitational forces are negligible. Comparison of modeling results
using the transit time equation, models for unsteady state flow, and
unsaturated flow suggest that the transit time equation underestimates
required liner thickness as a function of design life.
The proposed EPA transit time equation is based on steady state fully
saturated flow in the liner. GCA has modified this equation to include the
effects of reduced effective porosity, and negative pressure head (suction) on
the liner bottom. This modified equation still fails to capture the major
dynamics of the infiltration event of a wetting front moving into an
unsaturated soil. Primarily, this error is a result of ignoring capillary
forces during the transient infiltration period. This equation underestimates
design liner thickness required for containment.
The Green-Ampt (1911) model of infiltration has been suggested as a liner
design equation (see Pope-Reid Associates, 1982). This model approximates the
infiltration wetting front as a square wave to which saturated Darcy flow
analysis is applied. The resulting liner thickness equation captures more of
the dynamics of the infiltration event, but is difficult to apply because
unknown parameters must be specified, to which the solution is very
sensitive. In the presented comparison, this model led to both
over-conservative and insufficient design thickness.
The unsaturated flow governing equation can be linearized. Analytical
solutions are available for this linearized form, and liner thickness can be
iteratively determined. Although capturing much of the dynamics of the
wetting front movement, this technique is not recommended because liner
conductivity and moisture retention are assumed to be independent of moisture
content, and the constant parameters K* and D* must be determined by some sort
of fitting procedure, which cannot be performed a priori.
28
-------
RECOMMENDATIONS
Numerical simulation is recommended as a method for determining liner
behavior during infiltration. The nonlinear unsaturated flow governing
equation can be approximately solved, retaining the system nonlinearity, as
well as simulating both gravity and capillary forces. The finite element
method is particularly well suited to this application because it is more
accurate at representing the nonlinearities in hydraulic conductivity and
moisture retention across the discretized domain. The numerical model can
then be applied to particular liner simulations, and can simulate an endless
variety of soil types and moisture conditions. The numerical code is verified
(programming errors) by solving simple idealized problems for which solutions
are available. By changing input soil characteristics and boundary
conditions, the model is used to investigate and expose patterns of liner
behavior which can be incorporated into empirical relations of field test data
and long-term liner performance.
The recommended alternative liner thickness determination method requires
a combination of laboratory and field tests, and the use of an empirical
method for scaling test behavior to long-term liner reliability. Specific
tests should be conducted in the field and the lab on the undisturbed liner,
and its underlying layer. The results of these tests can be extrapolated to
facility scale behavior using relationships determined by physical and
numerical experimentation.
Numerical simulation techniques, such as finite difference and finite
elements, can be applied to the unsteady, unsaturated flow equation. These
solutions are not limited by the exclusion of complicated flow and porous
medium properties and can include a rigorous, albeit approximate,
representation of the system nonlinearity. Simulation of a wide range of
liner types and moisture conditions will lead to generic expressions, fitted
to numerical results, which will allow results of lab and field tests to be
scaled to the facility size and design life.
The approach demonstrated by Hamilton et al. (1981) provides an accurate
estimate of liner life, neither overestimating nor underestimating required
clay thickness. Laboratory tests of the type shown in Figure 12 of Section 3
or of the type shown in Figures 13 and 14 of Section 3 are appropriate. The
numerical model provides a means to scale up from the laboratory scale to
field conditions. Due care will be required in specifying details of
laboratory testing to duplicate essential features of field conditions such as
kneading action and permeant composition.
LINER SPECIFICATIONS
Guidance should be available to the prospective owner of storage
impoundment to specify practical aspects of selecting, and testing soils, and
installing and compacting liners. Based on our review we recommend use of the
following criteria as a minimum:
29
-------
Soil selection
- minimum plasticity index = 10 for clayey till, clayey fine
sand, or clay
sand lenses in borrow area to be less than 1 inch in thickness
Liner placement
- place liner in the presence of a qualified soils engineer
place and compact soil slightly wet of optimum (for heavy clays
optimum water content is slightly lower than the plastic limit)
achieve suction pressure
-------
SECTION 5
REFERENCES
Bear, J. Hydraulics of Groundwater. McGraw-Hill, New York, 1979.
Bruch, J. L., Jr. Finite Element Solutions for Unsteady and Unsaturated Flow
in Porous Media. California Water Resources Center, Contribution No. 151.
1975.
Brutsaert, W. Universal Constants for Scaling the Exponential Soil Water
Diffusivity, Water Resources Research, 1979, pp. 481-483.
Clapp, R. B., and G. M. Hornberger. Empirical Equations for Some Soil
Hydraulic Properties, Water Resources Research 14(4), 1978, pp. 601-604.
Crank, J. The Mathematics of Diffusion. Clarendon Press, Oxford, 1975, p.
389.
Croney, D., and J. D. Coleman. Soil Structure in Relation to Soil Suction
(pF). Road Research Laboratory, Department of Scientific and Industrial
Research. Published in the Journal of Soil Science, Vol. 5, No. 1, 1954. pp.
75-84.
Daniel, D. E., personal communication, August 25, 1982.
Day, Arthur. Correspondence to GCA from U.S. Environmental Protection Agency
Office of Solid Waste and Emergency Response. June 1982.
Elzeftawy, A., and K. Cartwright. Evaluating the Saturated and Unsaturated
Hydraulic Conductivity of Soils, in Permeability and Groundwater Contaminant
Transport ASTM STP 746, T. F. Zimmie and C. 0. Riggs, eds., American Society
for Testing and Materials, 1981, pp. 168-181.
Green, W. H. and G. A. Ampt. Studies in Soil Physics I: The Flow of Air and
Water Through Soils, Journal of Agricultural Science 4, 1911, pp. 1-24.
Gruber, P. A. Simplified Method for the Calculation of Unsaturated Hydraulic
Conductivity. Presented at AGU Spring Meeting, Philadelphia, PA,
May 31-June 4, 1982.
31
-------
Hamilton, J. M., D. E. Daniel, and R. E. Olson. "Measurement of Hydraulic
Conductivity of Partially Saturated Soils," Permeability and Groundwater
Contaminant Transport, ASTM STP 746. T. F. Zimmie and C. 0. Riggs, Eds.,
American Society for Testing and Materials, 1981, pp. 182-196.
Johnson, T. M., S. A. Rojstaczer, and K. Cartwright. Modeling of Moisture
Movement Through Layered Covers Designed to Limit Infiltration at Low-Level
Radioactive Waste Disposal Sites, Presented at AGU Spring Meeting,
Philadephia, PA. May 31-June 4, 1982.
Maslia, M. L., and R. H. Johnston. Simulation of Ground-Water Flow in the
Vicinity of Hyde Park Landfill, Niagara Falls, New York, U.S.G.S. Open-file
Report No. OF82-0159, April 1982.
Mclntyre, D. S., R. B. Cunningham, V. Vatanakul, and G. A. Stewart. Measuring
Hydraulic Conductivity in Clay Soils: Methods, Techniques, and Errors, Soil
Science 128(3), pp. 171-183, 1979.
McWhorter, D. B., and J. D. Nelson. Unsaturated Flow Beneath Tailings
Impoundments. J. Geotech. Eng. Div. ASCE GT(ll), 1979, pp. 1317-1334.
Miller, R. D., and E. Bresler. A Quick Method for Estimating Soil Water
Diffusivity Functions. Soil Science Soc. Am. J. 41, 1977, pp. 1020-1022.
Milly, P. C. D. Moisture and Heat Transport in Hysteretic, Inhomogeneous
Porous Media: A Matrix Head-Based Formulation and a Numerical Model, Water
Resources Research, 18(3), 1982, pp. 489-498.
Moore, Charles A. Landfill and Surface Impoundment Performance Evaluation
Manual. Submitted to the U.S. Environmental Protection Agency, Office of
Water and Waste Masnagement, by Geotechnics, Inc. SW-869, September 1980.
Ogata, A., and R. B. Banks. "A Solution of the Differential Equation of
Longitudinal Dispersion in Porous Media," U.S. Geol. Survey Prof. Paper No.
411-A, 1961.
Pinder, G. F., and W. G. Gray. Finite Element Simulation in Surface and
Subsurface Hydrology, Academic, New York, 1977.
Philip, J. R. The Theory of Infiltration, 6, Effect of Water Depth Over
Soil. Soil Science, 85(5), 1958, pp. 278-286.
Pope-Reid Associates, Inc. Correspondence to Mr. Arthur Day, U.S.
Environmental Protection Agency, February through April 1982. Incorporated
into GCA Work Assignment 68-02-3168, Task 72-(3), as Attachment B.
Trautwein, S. J., D. E. Daniel, and M. W. Cooper. A Case History Study of
Water Flow Through Unsaturated Soils. Presented at AGU Spring Meeting,
Philadelphia, PA, May 31-June 4, 1982.
32
-------
APPENDIX B
PARTIAL LIST OF AVAILABLE UNSATURATED FLOW MODELS
Brutsaert, W. F., A functional iteration technique for solving the Richards
equation applied to two-dimensional infiltration problems, Water Resource
Research 6(7), 1583-1596, 1971.
2-D INFILTRATION MODEL. Non-steady, vertical cross-section, up to 5 soil
types, user-oriented, FDM. 1971.
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.
Finite difference on 2-D unsaturated equations and application to well
hydraulics.
Cooley, Prediction of transient or steady state hydraulic head
distribution in unsaturated, anisotropic, heterogeneous, two-dimensional,
cross-sectional flow systems.
Corey, A. T.
AGDRG and WADSOR. Solves the two-dimensional, transient flow drainage
problems in unsaturated/saturated soils.
Davis, L. A. Computer analysis of seepage and groundwater response beneath
tailings impoundments, National Science Foundation, (available through
NTIS: NSF/RA-800054) Washington, D. C., 1980.
SEEPV. A transient flow model to simulate vertical seepage from a
tailings impoundment, including saturated/unsaturated modeling of
impoundment with liner, and underlying aquifer, 1980.
Dutt, G. R., M. J. Shaffer, and W. J. Moore, Computer simulation model of
dynamic bio-physicochemical processes in soils, Universtiy of Arizona
Agricultural Experiment Station Technical Bulleten 196,1972.
SALT TRANSPORT IN IRRIGATED SOILS. A transient one-dimensional, vertical
simulation of solute transport in the unsaturated zone, coupled with a
chemistry model, 1976.
B-l
-------
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, February 1969.
1-D finite difference model, over review of previous models.
Freeze, R. A., Three-dimensional, transient, saturated-unsaturated flow in a
groundwater basin, Water Resources Research, 7, 347-366, April 1971.
3-D finite difference and application to several flow problems.
Giesel, W., M. Renger, and 0. Strebel, Numerical treatment of the unsaturated
water flow equation: comparison of experimental and computed results,
Water Resources Research, 9, 174-177, February 1973.
Finite difference technique on pressure and experiment on sand column,
one dimensional.
Gillham, R. W., A. Klute, and D. F. Heermann, Hydraulic properties of a porous
medium: measurement and empirical representation, Soil Science Society
of America Journal, 40, 203-207, March-April 1976.
Presents an emperical extension to King's [1965] hysteretic curve fitting
model.
Haverkamp, R. and M. Vauchlin.
SIMTUS. 1-D, non-steady, unsaturated flow in isotropic soils, FDM. 1977.
Huyakorn, P.
SATURN2. Studies transient, two-dimensional variably saturated
flow and solute transport in anisotropic, heterogeneous porous media,
1982.
INTERA Environmental Consultants, Inc.
HYDROLOGIC CONTAMINANT TRANSPORT MODEL. A three-dimensional model to
simulate flow and solute transport in a saturated/unsaturated,
anisotropic, heterogeneous aquifer system, 1975.
Kaszeta, F. E., C. S. Simmons, and C. R. Cole
MMT-1D. Simulates transient, one-dimensional movement of radionuclides
a"nd other contaminants in saturated/unsaturated aquifer systems.
Khaleel, R. and D. L. Redell, Simulation of pollutant movement in groundwater
aquifers, Texas Water Resources Institute, Texas A&M University,
Technical Report No. 81, 1977.
B-2
-------
A two-dimensional vertical model for the simulation of unsteady two-phase
flow and dispersion in saturated-unsaturated porous media.
Konikow, L. F. and J. D. Bredehoeft, Computer model of two-dimensional solute
transport and dispersion in groundwater, in Techniques of Water Resource
Investigation, Book 7, Chapter C2, 1974.
Kraeger-Rovey, C. E.
LINKFLO. Simulates three-dimensional steady and unsteady saturated and
unsaturated flow in a stream-aquifer system.
Marino, M. A.
INFILTRATION FEM. Simulates transient movement and distribution of a
solute (introduced as a constituent or artificial recharge) in a
saturated-unsaturated porous medium.
McCracken, G.
MULTIPURPOSE. Solves any of the equations generally encountered in
subsurface flow and transport.
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.
Develops 1-D vertical finite element model and verifies by comparison to
infiltration solutions. Good review of current soil moisture research.
Molz, F. J.
2-D horizontal-saturated zone, 1-D vertical unsaturated zone, non-steady,
phreatic, finite-difference model, lumped, unsaturated zone. Auburn
University, 1974.
Narasimhan, T. N. and S. P. Neuman
FLUMP. 2-D in vertical or horizontal plane, phrestic, confined and
leaky, specific storativity, versatile, flexible, tested, finite-element
model, U. Berkeley, 1975.
National Energy Software Center (NESC)
ODMOD. Prediction of coupled one-dimensional, vertical movement of water
and trace contaminants through layered, unsaturated soils.
B-3
-------
Neuman, S. P., Saturated-unsaturated seepage by finite elements, Journal of the
Hydraulics Division, ASCE, 99(HY12), 2233-2250, December 1973.
UNSAT II. Computes hydraulic heads, pressure heads, water content,
boundary fluxes and internal sinks and sources in a
saturated/unsaturated, nonuniform, anisotropic, porous medium under
nonsteady state conditions.
Pickens, J. F.
UNFLOW. Simulation of two-dimensional (cross-sectional) transient
movement of water in saturated-unsaturated nonuniform porous media, 1977.
Pikul, M. F., R. L. Street, and I. Remson, A numerical model based on coupled
one-dimensional Richards and Boussinesq equations, Water Resources
Research, 295-302, April 1974.
Quasi 2-D finite difference model: horizontal flow in saturated zone,
vertical flow in unsaturated zone.
Reed, J. E.
2-D horizontal-saturated zone, 1-D vertical-unsatui .it^d one, non-s>ceady
phreatic, confined and leaky, specific storativity, roots, river,
aquitard, user oriented, mass storage, finite-difference model. USGS,
1974.
Reeves, M. and J. 0. Duguid, Water movement through saturated-unsaturated
porous media: a finite-element galerkin model, Oak Ridge National
Laboratories, Report No. ORNL-4927, February 1975.
User's manual and development of 2-D vertical flow model, comparison to
Freeze's (1971) simulations.
MOISTURE TRANSPORT CODE. A two-dimensional transient model for flow
through saturated/unsaturated porous media.
Reisenauer, A. E.
TRUST. Computes steady and nonsteady pressure head distributions in
multidimensional, heterogeneous, variably saturated, deformable porous
media with complex geometry.
Robertson, J. B.
TRA POND MODEL. Simulates subsurface transport of radionuclide solutes
from seepage pond through perched water zones to regional aquifer.
Selim, H. M., R. S. Mansell, and A. Elzeftawy, Distributions of 2,4-D and
water in soil during infiltration and redistribution, Soil Science,
121(3), 176-183, 1976.
B-4
-------
NMODEL. Predicts water and nitrogen transport and transformations under
transient and steady unsaturated water flow in homogeneous or
multilayered soils.
Simmons, C. S., see: User's Manual Unsaturated groundwater flow model -
UNSAT1D, EPRI Report CS-2434-CCM.
UNSAT1D. One dimensional simulation of unsteady vertical unsaturated
flow, 1978,
Skaggs, R. W., Combinationsurface-subsurface drainage systems for humid
regions, J. of the Irrigation and Drainage Div., ASCE, 106(lR4), 265-283,
1980.
DRAINMOD. An unsteady, one-dimensional, horizontal/vertical,
saturated/unsaturated model to simulate watertable position and soil
water regime above water tatile for artificially drained soils.
Terry, J. E.
SUPERMOCK. Simulates transient stress and response in a
saturated-unsaturated ground water flow system including a water-table
aquifer overlying a confined aquifer.
Washburn, J. F.
MMI-DPRW. Predicts the transient three-dimensional movement
of radionuclides and other contaminants in unsaturated/saturated aquifer
systems.
Wierenga, P. J., Solute distribution profiles computed with steady-state
and transient water movement models, Soil Science Society of America
Journal, 41, 1050-1055, 1977.
Compare predicted solute concentrations using steady vs transient
unsaturated flow. Under pulsed boundary conditions, both flow models
produced very similar concentration histograms. Silty clay loam.
Exponential moisture characteristic and relative permeability.
Wind, G. P. and W. Van Doorne, A numerical model for the simulation of
unsaturated vertical flow of moisture in soils, Journal of Hydrology, 24,
1-20, 1975.
Yeh, G. T. and D. S. Ward, FEMWATER: a finite-element model of water flow
through saturated-unsaturated porous media, Oak Ridge National
Laboratories, Report No. ORNL-5567, October 1980.
User's manual and modifications to 2-D vertical flow model of Reeves and
Duguid (1975). Based on pressure head.
B-5
-------
APPENDIX C
LISTING OF EXAMPLE COMPUTER PROGRAM
SOILINER
C-l
-------
09/26/83
08:39:59
GCAGD.SOILINER.FORT.DAT*
uuuu an u
00000020
OOOOOC30
00000010
00000050
OOOOOC60
00000070
00003080
00000090
00000100
00000110
30000120
n fi n n f 1 1 -in
UUUUUIJU
00000140
00000150
00000160
00000170
00000180
00000 190
00000200
30000210
00000220
00000230
00000210
00000250
00000260
00000270
00000280
00000290
00000291
00000300
00000310
00000320
00000333
00000310
00000350
00000360
00000370
30000380
00000390
00000100
00000110
0000012C
00000130
00000110
00000150
00000160
nnnnni»"7fi
UUUUUHfU
00000180
nnrirt n A Q n
UUUUUH^U
00000500
00000510
00000520
00000530
00000510
00000550
00000560
00000570
COOOO 580
00000590
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
CBUG
C
C
C
C
C
C
C
r
l*
C
C
SOILINER. FORTRAN FINITE DIFFERENCE OF VERTICAL UNSATURATED
INFILTRATION
DAN GOODE JUNE 1983
COPYRIGHT 1983 GO/TECHNOLOGY DIVISION
213 BURLINGTON ROAD
DEDFORD, MH 31730
(617) 275-5111
COMMON/MASSb/STOH.TOTVl ,TOTV2,FLUX10iFLUX20
COMMON/DEVICE/ IRQ »IPRT, IFILE
COMMON/ FILES /I FGRDt IFSUIL t IFINIT . IFPOUT. IFMOUT* IFSOUT . IFLUX » IFZOUT
COMMON/ INF 0/NUMNP .NUMEL tNHMl « ,^PM2t NUhEL2 » AM It IbP ACE .1 TOL.
* DEPTH,ENDTIM,OTfDTHAX»ALPHA»PSIBOT.H.NZOUT«
* PSIMIT,NCPTS,ERKMAXfMAXIT.CiHPARMtMAXNTt ISTEOY
COMMON/TIMES /TIME tTIMEl .NT« ER RtNOUTiTOUT < 1 0 > t NOUT1, TOUT1
COMMON/ ERRORS/ 1 ERR. I WARN
DIMENSION PSK200).PSIOLD(200).02(200).C(200).
* RM200) »*ORH<200) «RMLURV(1C> t
* RKCURVC i:>,FSICRV(ia>tRHOISTl200)t
* STIFF(3,200),DPSI<203) ,F(200)t
* Z(200>t SATK»200>«PSICRT<200> .POR4200)
DIMENSION VNEU (200) tVCL 0(200) . OLDC (20 0 ) »DZN<200> .STARK (200) t
* IZOUT(lC),ISOIL(200)t3ETA(200)
DIMENSION RKL(100)tCL(100)«RKSTL(130),OLOML(100)
DATA NMAX/2TO/
IUARN=0
IRD = 5
IPRT=7
READ AND WRITE PROBLEM PARAMETERS AND INITIALIZE
CALL INPUT
-------
00000600 C
00000610 C
00000620 120
00000630 C
00000640 C
00000650
ITER=ITER+1
COMPUTE MOISTURE CONTENT AND CONDUCTIVITY
CALL SPROPl (PSI»PSICRV,RMCURV ,*KCURV,RMOIST,C»
00000660 * RK,R^STL,CL,RKL, STARK, SA TK.POR . I SOIL .PS1CRT,
00000670 C
oooooteo c
00000690
00000700 C
00000710 C
00000720
00000730 C
00000740 C
00000750
00000760 C
00000770 C
00000780
00000790 C
00000800 C
00000810 C
00000820
00000830 C
00000840 C
00000850
00000860
00000870
00000880 125
00000890 C
00000900 CBU6
00000910 C
00000920 C
00000930 C
00000940
00000950
00000960
00000970 C
nnnnficfin f
UuuU U 70 U L
00000590 C
nnnf*ir*nn f
UUUUluUU \,
00001010 C
00001C20
00001030
00001040 C
00001050 C
00001060 C
00001070 140
00001080
00001090 C
00001100
00001110
00001120 C
flflrtfiii'^fi (*
U U UU 1 4. J U t,
00001140 C
n ft n n 1 1 c* n c
U U UU 1 ID U C
00001160 C
00001170 C
00001180
00001190
00001200 C
00001210
ZERO STORAGE VECTOR
CALL CLEAR(C.NUKNP)
BUILD VELOCITY VECTOR
CALL BUILD VC VNEW.PSI ,STARK,DZN,DZ,BET A ,RK )
BUILO STIFFNESS MATRIX
CALL BUILDS(STIFF,RK,STARK,C,OLDC,DZN,DZ,BETA)
BUILD FORCING VECTOR BOUNDARY CONDITION TERMS
CALL BUILDF(F,C,OLDC»PSI,PSIOLO, VNEU.VOLD)
SOLVE MATRIX EQUATION USING THJMAS ALGORITHM
RESULT IS INCREMENTAL CHANGE IN PSI, OPSI
C«LL THOMAS(STIFF»F,DPSI,WQRK,fJPM2,0)
UPDATE PSI
DO 125 I=1,NPP2
N = I*1
PSI GOTO ItO
IFUTER.LE.MAXIT) SOTO 120
END OF MAIU ITERATION LOOP
URITE (IPRT.2000) NT
IWARN=IUARN+1
URITE REQUESTED OUTPUT TO PRINT AND PLOT FILES
CALL SPROPl (PSI,PSICRV,RMCJRV ,RKCURV,KMOIST,C,
* RK,RMSTL,CL,RKL,STARK,SATK,POR,ISOIL,PSICRT,
ALPHA=ALPHAS
AM1 = 1. DO-ALPHA
END OF STEADY STATE SOLUTION
OUTPUT INITIAL CONDITIONS AND/OR STEADY STATE. SOLUTION
1)
CALL OUTPUT! PSI, RMO I ST,RK,Z,RKL,CL,RM STL, STARK,POR,DZ, ITER)
IFCNZOUT.GT. 0) CALL ZOUT< IZOUT»PSI ,RM OIST )
IF(MAXNT.LE.O) GOTO 50
C-3
-------
00001220 C
00001230 C
0000121C
00001250 C
C0001260 C
30001270
00001280 C
00001290 C
nnrtni"?rn p
UOOOlJUu L
00001310 C
o n n rt 1 "ion f
u U UU 1 -c u U
00001330 C
00001240
00001350 10
00001360
00001370 C
00001380 C
00001290
00001*00 *
00001110 C
00001420 C8UQ
00001430 C
00001440 C
00001450 C
00001460 C
00001470 C
00001480
00001490 C
00001500 20
00001510 C
00001520 C
00001530
00001540 •
00001550 C
00001560 C
00001570
00001580 C
00001590 C
00001592 CBUG
00001593 6&63
00001600
C0001610 C
00001620 C
00001630
00001640 C
00001650 C
00001660 C
00001670
00001660 C
00001690 C
00001700
00001710
00001720
00001730 25
00001740 C
00001750 CBUG
00001760 C
00001770 C
00001780 C
00001790
00001800
00001810
BUILD VELOCITY VECTOR IF NOT ALREADY COMPUTED FROM STEADY STATt
IF(ISTEDY.NE.l) CALL PUILDV ( V NE J.PSIf S TARKt DZN.JZtB ET A tRK )
READ IN NEW BOUNDARY CONDITIONS IF INITIAL CONDITION IS SS
IFilSTEDY.tG.l) CALL NEyBC«PSI>
START TRANSIENT aOLUTIuN
CALL MASBALCPSI.RMSTL.OLDML .STARKt DZ.-l)
ITER=0
NT=NT+1
SAVE OLD VECTORS AND PREDICT NEW PSI FROM V
CALL PREDCT(PSI»PSIOLD,VNEWtVOLD,C.OLDC?
OLDML.RMSTL)
CALL OUTPUTtPSI,R«OIST,RK,Z»RKL«CL.RMSTL»STARKtPOR.DZtITER>
MAIN ITERATION LOOP
TIP!E = TIME1+DT
ITER=ITER+1
COMPUTE MOISTURE CONTENT AND CONDUCTIVITY
CALL SPROPKPSI tP SICR Vf RMCURV ,3KCURV .RMOIST.C t
RK,fil*STL,CL.RKLt bTARK,SATK,POR t ISOIL.PSICRT, 0)
BUILD VELOCITY VECTOR
CALL BUILOVf VNEU, PSI « STARK ,DZNtOZ t BETA »RK>
BUILD STIFFNtSS MATRIX
URITE(b«6663) NT.OT.DTMAX .T IME.ENDTIM
FORMAT( I5.1P4E15.4)
CALL BUILDS(STIFF,RKtSTARKt Ct OLDCtDZN.DZ tBETA )
BUILD FORCING VECTOR TRANSIENT AND BOUNDARY CONDITION TERMS
CALL BUILDF(F,C,OLOC.PSItPSIOLO,VNEU,VOLD)
SOLVE MATRIX EQUATION USING THOMAS ALGORITHM
RESULT IS INCREMENTAL CrIAM&E IN PSIt QPSI
CALL THOMASC STIFF tf t DPS I, yORK «NPM2 « 0)
UPDATE PSI
DO 25 I=1.NPM2
N=I + 1
PSI(N)=PSI(N1»OPSI(I>
CONTINUE
CALL OUTPUT
-------
00001820
n n n n 1 Q i n
U U U U 1 oj U
00001(340
00001 850
00001860
0000187C
00001880
G0001b90
00001900
00001910
OC00192C
00001930
00001940
00001950
00001960
0000197C
00001 S80
00001990
00002000
00002010
00002C20
00002030
00002040
00002050
00002060
00002070
00002C80
00002090
00002100
00002110
G0002120
C0002130
00002140
00002160
n n nn o i 7 n
U U UU c 1 * U
00002180
00002190
000022CO
OC002210
00002220
00002230
00002240
00002250
00002260
00002270
00002280
00002290
00002300
0000231P
00002320
00002330
00002340
00002350
00002360
00002370
00002380
00002390
00002400
00002410
30002420
00002430
C
C
C
C
C
C
C
30
35
C
C
C
C
C
C
C
50
C
2000
2010
C
C
C
C
C* • • •
C
C
C
C
END OF MAIN ITERATION LOOP
WRITEC IPRT.2000) NT
IWARN=IWARN+ 1
WRITE REQUESTED OUTPUT TO PRINT AND PLOT FILES
CALL SPROPKPSI,PSICRV,RMCJRV,RKCURV,KMOIST,L,
* PK,RMSTL,CL»RKL»STARK,SATK,POR,ISOIL«PSICRT, 1)
IF(NZOUT.GT.O) CALL ZOUT< IZOJT.PSI »RMOIST )
IF
IFCNOUT1.GT.NOUT.OR.TGUT1 .LE. O.DO) NOUT=0
CONTINUE
IF(UOUT.EQ.-l.OR.TIME.GE.ENDTIM.OR.NT.GE.MAXNT)
* IOUT=1
IF(IOUT.EQ.l) CALL OUTPUT (PS I ,RMOIST ,RK,Z ,RKL ,CL ,RMSTL ,STARK,POR,
* D2.ITER)
CALL MASBAL(PSI,RMSTL,OLDML,!>TARK,DZ,IOJT)
IFITIME.&E.ENDTIM.OR.NT.GE.MAXNT) GOTO 50
NEW TIMES
TIME1=TIME
IOUT=0
CALL NEWTIM«PSI,PSIOLD, IOUT)
GUTO 10
END OF SOLUTION PERIOD
WRITE(IPRT,2010) IJARN
STOP
FORMATt/* ***» WARNING ***• MAXIMUM ITERATIONS EXCEEDED *,
*»TIME STEP 4,I5)
FORMAT
LOOP OVER ELEMENTS AND COMPUTE SOIL PROPERTIES AT NODES
AT EACH END OF ELEMENT. SINCE ADJACENT ELEMENTS ,1AY HAVE
DIFFERENT SOIL TYPES, A NODE HAS TWO VALUES FOR EACH PROP.,
ONE FOR ELEMENT ON TOP, AND ONE FOR BOTTOM ELEMENT.
COMMON/MASSB/STOR,TOTV1,TOTV2,FLUX10,FLUX20
COMMON/DE»/ICE/IRD,IPRT,IFILE
COMMON/INFO/NUMN?,NUMEL,NP«1,NPM2,NUMEL2,AM1,ISPACE,ITOL,
* DEPTH,ENDTIM,DT,DTMAX,ALPHA,PSIBCT,H,NZOUT,
* PSINIT,NCPTS,ERRMAX,MAXIT,CHPARM»MAXNT,ISTEDY
C-5
-------
00002440
00002450
00002460
00002470
00002480
00002490
00002500
00002510
00002520
00002530
00002540
00002550
00002560
00002570
0000258C
00002590
00002600
00002610
00002620
00002630
00002b40
00002650
00002660
00002670
00002680
00002690
00002700
00002710
00002720
00002730
00002740
00002750
00002760
00002770
00002780
00002790
00002BOO
00002810
00002820
00002630
00002840
00002850
00002860
00002870
0000288C
00002690
00002900
00002S10
00002520
00002930
00002940
00002950
00002560
00002570
00002580
00002990
00003000
00003010
00003020
00003C30
00003040
000,03050
C
C
C
C
CBUG
6666
C
8
CBUG
C
C
C
22
CDP
C
C
C
C
C
C
C
C
C
C
COMMON/ERRORS/IERR,I WARN
COKMON/FILES/IFGRD,IFSOIL,IFINIT,IFPOUT,IFMOUT,IFSOUT,IFLUX.IFZOUT
DIMENSION PSI(l>tPSICRV(l),RMCURV(l)f
RNOISTC1),RKL<1>,SATK(1),PuR(l), CU> ,
RKCUhV(l),RK(1),STMRK(1),CL<1>,RMSTL<1>,I SOIL<1),PSICRT(1)
LOOP AND CALCULATE SOIL PROPERTIES
00 10 L=1,NUMCL
USOIL=ISOIL(L)
TOP NODE
IP=IP+1
URITE<6»6666> L»K,N,IP,JSOIL,PSI(N),PSICRT(L)
FORM»T(5I5,2F10.2)
IF(PSKN).GT.PSICRT (L » GOTO 5
GOTO 22
BOTTOM NODL
IP=IP+1
K-2
URITE(6,6666> L,K,N,IP,JSOIL.PSI(N),PSICRT(L)
IF(PSICN).GT.PSICRT
RLTUHN
END
SUBROUTINE BUILOV (V, P SI .STARK ,DZN, DZ, BETA ,RK )
COMPUTE C*D
-------
C0003680
30003690
0000370C
00003710
P0003720
30003730
00003740
D000375C
00003760
00003773
0000378C
00003790
00003800
C0003610
00003820
30003830
00003840
00003650
00003860
00003870
00003880
00003890
00003900
00003910
00003920
00003930
00003940
00003950
00003960
00003970
00003980
00003990
00004000
00004010
00004020
00004030
00004040
00004050
00004060
00004070
00004080
00004090
00004100
00004110
00004120
00004130
00004140
00004150
00004160
00004170
00004180
00004190
00004200
00004210
00004220
00004230
0000424C
00004250
0000426C
00004270
00004280
00004290
C THIS VERSION HAS SPECIAL ALGORITHM FUR VARIABLE SPACING
C
CCMMON/MASS6/STORtTOTVl»TOTV2»FLUX10tFLU*.>0
COMMON/TIMES /TIME »T I «tl,NT» ERR »NGUTtTOdT( 1C) ,
COMMON/DEVICE/ IRO,IPRT, IF ILd
COMMON/ I NFO/NUMNP ,MUM CL ,NF-« 1 , NPM2 , NUMEL2 , AMI ,
NOUT1,TOUT1
ISPACE,ITOL,
* DEPTH, ENDTIM,OT,OTM AX, ALPHA, PSItOT,H,NZ OUT,
* PSINIT,NCPTS,£RR*AX,MAxIT ,Crt=>AF, PS I«l>, STARK (1),DZN(1),9Z(D»RK(1>»BLT All)
C
IF(ISPACE.EG.l) GOTO 20
C
DO 10 I=1,NPM2
L = I
LP1=L*1
N - I + 1
C
V(I) = ( STARK(LPl) * (PSI(N+1)-PSI(N>)/ DZ(LP1
* STARML) * ( PSI < N)-PSI (N-l) )/ DZ(L) »
* STARMLP1 >-STARK(LM /DZN(I»
10 CONTINUE
GOTO 50
C
C VARIABLE SPACING ALGORITHM
20 DO 40 I=1,NPM2
L-I
LP1=L+1
N = I + 1
NH1=H-1
NP1=N+1
B1=BETA(I)
IF(B1 .EQ.l.DO) GOTO 30
B2=1.UO/B1
DZNIrU2N(I)
VtI) = iB2*STARK«LPl)*«PSHNPl)-PSI(N»/DZ(LPa)
* bl*STARK(L)*(PSHN>-PSI (NMD )/DZ(L) «
* (B1-B2) *RK(N>*( (PSI(W>1> -PSKNMl ) )/(DZNI
* B1*STARK(L) +D2*iTARK(LPl) ) /JZMI
GOTO 40
) -
-
+DZNI) +1) -
30 V(I) = ( STARKCLP1) * ( PS II NP 11 -PS I( N)> / DZ(LPl) -
* STARK(L) * (PSKNI-PS1 (NMD)/ DZ(L) +
* STARKCLP1 )-STARK(L) ) /DZN(I>
4C CONTINUE
C
50 CONTINUE
CBUG CALL OUMPK'V BUI LDV • ,V ,NPM2)
C
999 RETURN
END
C
SUBROUTINE PREDCT(PSI ,PSI OLD. VNE 1 , VOLD ,C , OLDC
* DLOML,RMSTL >
C
C
C.... STORE OLD VECTORS AND PREDICT NEW PSI
C
COMM3N/MASSB/STOR ,TOTV1,TOTV2»FLUX10,FLUX20
CGKMON/DEV1CE/IRD.IPRT,IFILE
COMMON/ I NFO/NUMNP,NUMLL,NPM1,NPM2,NUMEL2, AMI,
,
ISPACE,ITOL,
C-8
-------
00004300
00004310
C0004320
0000*330
00004340
00004350
00004360
00004370
00004380
00004390
00004400
00004410
0000442C
00004430
COOC4440
00004450
00004460
00004470
30004160
00004490
0000450C
p.nrtftACl n
UUUUHO1U
D0004520
00004530
00004540
nfinni^cn
UUUUH ^30
00004560
00004570
00004580
00004590
00004600
00004610
00004620
00004630
00004640
OC004fa50
00004660
00004670
00004680
30004690
00004700
00004710
00004720
00004730
00004740
00004750
00004760
00004770
00004780
00004790
00004800
30004810
00004&20
00004630
30004840
00004850
00004860
00004870
00004880
00004690
00004900
00004910
C
C
10
C
CBUG
C
C
C
C* • • •
C
C
C
C
C
C
C
C
10
C
C
C
20
* DEPTH,ENOTIM,DT,DTMAX,ALPHA,PSIBUT,H,r,ZUUT,
* PSINIT,NCPTS,ERMAX,MAXIT,CH:>AKM,MAXNT,ISTLDY
COMMON/TIMtS/TIMC.TIKEl .NT. LR K, NOUT ,T UUT ( 10) , NGUTl.TuUTl
COMMON/ ERRORS/ I ERR, I WARN
DIMENSION PS I(l>.PbIOLD(l)«VNEU( l),VOLu(l).C( 1) tOLDC( 1).
* OLDML(1),RM3TL(1)
CALL COPY(PSI.NUMNP.PSIOLD)
CALL C JPY( VNEU.NPM2. VCLO)
CALL COPYIC.NUMNP.OLDO
CALL COPY(RMSTL,NUMEL2,OLDML>
PREDICT NEU PSI FRO" V
DO 10 I=1,NPM2
N = I + 1
IF(OLDCCN) .LE.O.DO) GOTO 1 Cl
PSI/DT
STIFFl3.I)=-ADZINV»STARKtLPll /DZ(LPl)
CONTINUE
GOTO 50
VARIABLE SPACING ALGORITHM
DO 40 I=1.NPM2
L=I
C-9
-------
00001920
LP1=L+1
00001910
00001550
0000 1960
C0001970 C
00001980
00001^90
00005COC
00005C1C
00005020
G0005030
00005010
30005050
00005060
00005070 C
00005080 C
00005090 C
00005100 30
00005110
00005120
00005130
00005110 1C
00005150 50
00005160 CBUG
00005170 999
00005180
00005190 C
00005210 C
B0005220
00005230 C
00005250 C
00005260 C....
00005270 C
00005280
00005290
00005300
00005310
00005320
00005330
00005310
00005350
00005360 C
00005370
00005380
00005390
00005100
00-005110
00005120 10
00005130 CBUG
00005110
00005150
00005160 C
00005160 C
00005190
00005500 C
00005E20 C
00005530 C....
Bl^BETAU)
ADZINV=ALPHA/DZN< I)
IF(B1 .EQ.l.DO) GOTO 30
B2=1.DO/B1
N i*l 1 - N - 1
NPl^N+1
DZN2 = DZ (D+DZCLP1)
STIFfrU.I> = -ADZINV*/DZ«<32-Bl)«R.K(N)/DZN2>
STIFF(2.I)=ADZIfJV*(Bl*STARK(L)/DZ(L)»B2»STARK(LPl )/D2(LPl )) +
• « ALPHA*C«N)*A11*OLDC(N> )/UT
STIFF(3,I)--ADZINV* ( b2*ST AR K< LP 1 ) /bZ < LPI > + *RK( N
GOTO IP
CONSTANT SPACING FOR THIS NODE
SriFFUtI)--ADZINV*ST,JRK(L>/DZ
S7IFF(2,I)-ADZINV*(ST1RK(L)/3Z(L)+STARK(LP1)/OZ(LP1» +
» ( ALPHA*C(N)»A11«OLDC(N) )/iT
STIFF(3»I)--ADZINV*STARK( LP 1 ) /OZ (LH)
CONTINUE
CONTINUE
CALL DUMPK'S BUI LOS» ,STI FF ,NPM2«3)
RETURN
END
SUBROUTINE bUILDF
BUILO FORCING VECTOR WITH BOUNDARY CONDITIONS AND TRANSIENTS
COMMON/MASSB/STOR,TOTV1,TOTV2.FLUX10,FLUX20
COMMON/DEVICE/IRD,IPRT,IFILE
COMMON/ I NFO/NUMNP .NUMEL «NPM 1 , NPM2 .NUMEL2 » AMI » ISP ACE .ITOL.
* DEPTH,ENDTIM,DTtDT"AX»*.LPHA,PSIBOT.HtNZUUT
»
* PSINITtNCPTS, ERRMAX.MAXITtCHPARM.MAXNT, ISTECY
COMMON/ERRORS/IERR.IJARN
DIMENSION F (1).CtVOLC(l>tOLDC(l)
DO 10 I=1.NPM2
N=I»1
F -
* +
* AM1*OLDC(N) )/DT
CONTINUE
CALL DUMPK'F BUILDF • ,F ,NP«( 2)
RETURN
END
SUBROUTINE NEUTI1 (KNEW, OLD, KODE )
COMPUTE NE» TIME STEP BASED ON CHANGE DURING PREVIOUS
STEP
C-10
-------
00005540
00005550
C0005560
00005S70
00005580
00005590
000056CO
00005610
0000562C
00005630
00005640
00005650
00005660
00005662
0000'5670
00005680
00005690
00005700
CC00571C
00005720
00005730
00005710
00005750
00005760
00005770
00005780
00005790
00005800
00005610
00005820
00005830
00005840
00005850
00005860
G0005870
00005b60
00005890
00005900
00005910
0 0 00 5 52 0
00005930
00005S40
00005950
00005960
00005970
00005980
00005990
00006000
00006C10
0000602C
00006030
00006C4C
00006050
00006C60
00006C7Q
00006080
00006090
00006100
00006110
00006120
00006130
00006140
C
C
C
C
CDP
10
C
C
C
C
20
C
C
C
C
~c
C* • • •
C
C
C
CDAN
COMMON/DF.VICE/IRD,IPRT, IF1LE
COMMON/TIMES/TIME.TIMLl,NT,ERR,NOUT,TOUT< 1 0) , NOUT1, TOUT1
COMMON/ I NFO/NUMNP ,NUM EL »NPM 1, NPM2, NUMEL2 , AMI , ISPACL ,1 TOL,
• DEPTH, END TIM,QT ,3T«AX,ALPHA,?S1BOT»H«NZOUT»
* PSINIT.NCF'TS.t^RlMX.MAXlT.CH'AR^MAXNT.ISTL'OY
COMMON/ ERRORS/ I ERR, I WARN
DIMENSION RNEu(l),OLD(l >
FIN3 LARGEST ABSOLUTE CHANGE
CHMAX=O.DO
DO 10 I=2,NUM£L
CHANGE NEXT CARD FOR DOUBLE/SINGLE PRECISION
CHANGE = ABS(RNEU(I)-OLD( I) )
IF,BET A( l),DZNtl),Z(l),?SICRT(l),
* IZOUT(l)
COMMON/MASSB/STOR,TOTV1,TOTV2,FLUX10,FLUX20
COMMON/ TIME S/TI ME, T I MTl, NT, ERR, NOUT, TOUT ( 1 0) , NOUT1, TOUT1
COMMON/ DE VIC E/IRD,IPRT, IF 1LE
COMMON/INFO/NUMNP,NUMEL,NPM1,NPM2,NUMEL2,AM1,ISPACE,ITOL,
* DEPTH,ENDTIM,QT,DT,1AX .ALPHA, PS IBOT,H ,NZ OUT,
PSINIT,NCPTS,ERRMAX,MAXIT,CHPARM,MAXNT,ISTEDY
COMMON /ERR OR S/IERR, I«ARN
COMMDN/FlLES/IFGRD,IFSOIL,IFINIT,IFPOUT, IFMOUT, IFSUUT , IFLUX , IFZOUT
TIME=O.DO
TIME1=O.DO
NT = 0
USE ITOL=0
ITOL=0
C-ll
-------
00006150 C
00006160
00006170
00006180 C
00006190
0000620C
00006210
00006220
00006230
00006240
00006250
00006260
00006270 C
000062BO C
00006290 C
OC006300
00006310
OC006320
00006330
00006310
00006350
00006260
00006370
0000638Q C
00006390 C
00006100 C
00006110
00006120
00006130
00006110 C
0000615C C
00006160
00006*70
00006180
00006190
00006500
00006510
00006520
00006530
00006510 C
00006550 C
00006560
00006570
00006580
00006590
00006600
00006610
00006620 C
00006630 C
00006610 C
00006650
00006660
00006670
00006680
00006690
00006700
00006710 C
00006120 C
00006730 C
00006710
00006750 C
00006760
READURD»1001> TITLE
URITEUPRT.20C1) TITLE
READ(IRD,1002) IFGRD»IFSOIL.IFINIT»IF1LE,IF POUT,IFMOUT«IFSOUT»
* IFLUX.IFZOUT
IF(IFGRD.LE.O) IFoROrlRQ
IFdFitML.LE.O) IFbOIL=IRD
II-UrINIT.LE.O> IFINIT-IRO
IF(IrILL.Lt:.0) IFILE = 10
URITEUPRT,20D2) IFvJRD, IF SOIL ,IF I NI T, IF ILE , IF POUT, I FMOUT. IF SOUT,
* IFLUX.IFZOUT
TEM3QRAL CONTROL PARAMETERS
READtIRD»1C03> IN ITAL ,1STEOY, tlNDTlM.DTtQTMAX,MAX NT,ALPHA,ERRMAX,
* MAXIT.CHPARH
IF(DT.LE.O.DO) DT=1.00
IF(£RHMAX.Lt .0.00) ERRMAX=Q.0500
IF(MAXIT.Lt.O) MAXIT=3
AH1=1.DO-ALPHA
yRITE NUMNP»NUMEL,IREGLR.ISPACE
REA3 SPATIAL OUTPUT PARH3
READIIRD,1030) NZOUT
URITE(IPRT,2006) NZOUT
IF(NZOUT.LE.0) GOTO 96
RLAOCIRD,10Cfc> (IZOOT(IZ) ,IZ = 1,MZOUT>
aRITE(IPRT,2007>
-------
00006770
00006780
00006790
JC006800 C
OC006G10
00006020
00006822
30006830
00006840 C
00006650
00006660
00006870
oooo&teo
C0006890
C000690C
00006510
00006920
00006930
00006940
00006550 C
00006960 C
00006970 C
00006980
00006990
00007000
00007010 C
00007020
00007030
00007C40
00007050
00007060
00007070
00007080
00007090
00007100
00007110 C
00007120 C
00007130 C
00007140
00007150
00007160
00007170
00007180 C
00007190
0000720C
00007210
00007220
00007230
00007240 C
00007250 C
00007260 C
C00072/0 C
00007280
00007290
00007300
00007310
00007320
00007330 C
00007340 C
00007350
00007360 C
00007370
READ IN,Z
IN1=IN2+1
IN2=IN1*NLZ-1
DO 14 IN=IN1,IN2
ZZ=ZZ-DZZ
ZtIN)=ZZ
14 CONTINUE
IFUN2.LT.NUMNP) GOTO 13
COM'UTE S'ATIAL DIFFERENCE TERMS
OVER ELEMENTS
47 DO 48 L = 1«NUME.L
DZ=ZlL»l>-ZtL>
48 CONTINUE
OVER INTERIOR NODES ONLY
DO 49 I=1.NPM2
N = I + 1
NH1=N-1
NP1=N+1
DZMI )-(Z I,DZ(I)
69 CONTINUE
WRITE*IPRT,2069) NUMNP,Z
DO 68 I=1,NUME.L
READ!IF SOIL,10701 IN,IS 01L(I),SATK(I),POR(I),PS1CRT (I )
JRITEtIPRT,2070) IN,I SOIL(I),SATK PSI.MIT
C-13
-------
00007380
00007390
00007400
00007110
00007120
00007130
C0007110
00007150
00007160
00007170
00007180
00007150
00007500
00007510
00007520
00007530
00007110
U0007550
00007560
00007570
00007580
00007590
00007600
00007610
00007620
00007630
00007(10
00007650
00007660
00007670
00007680
00007690
00007700
00007710
00007720
00007730
00007710
00007750
00007760
00007770
00007780
00007790
00007800
00007810
00007820
00007830
00007810
00007650
00007660
00007870
00007880
00007890
00007^00
00007910
00007920
00007930
00007910
00007950
00007960
00007970
00007980
00007990
URITE(IPRT.2C60> PSINIT
DO 52 I=1,NUMNP
PSKI J-PSINI T
52 CONTINUE
GOTO 55
C VARIABLE INITIAL CONDITION
53 DO 51 1=1,NUMUP
R£AD( IFINIT.1061) IN, PSKI)
51 CONTINUE
URIT£(IPRT,2090> (PS I (I >» I=1»NUHNP)
r-
U
C SET PRESSURE BOUNDARY CONDITIONS
C
55 CALL NEWBC(PSI)
r
999 RLTURIM
C
10C1 FORMATU8A4)
1002 FORMATOI5)
10G3 FORMAT<2I5t3F10.0,I5,2F10.0,I5,F10.0>
1006 FORMATU6I5)
1010 FORMAT(3I5>
1030 FORMAT( 15)
1010 FORMATC3F1C. 0)
1050 FoKMAT<8F10.0>
1061 FORM»T< 15, Fl 0.0)
1070 FORMATC2I5t3F10.0>
1080 FORMAT(I10,F 10.0)
1081 FORMAT(FIO.C)
C
2000 FORMAT!//1 TEMPORAL DISCRETIZATION PARAMETERS'//
* INITIAL CONDITION PAR M' « 1 DC 1H. ) , • Ih ITAL = ' , 1 1 O/
* IF INITAL EQ 1, READ INITIAL PSI FOR EACH NODE*/
« OTHERWISE, RExD CONSTANT INITIAL PSI'/
• STEADY STATE FARM ' , 12 ( 1 H. ) , MSTL& Y = • , 1 1 3/
* IF ISTEDY EQ 1, C01PUTE STCADY STATE SOLUTION'/
* OTHERUISE, COMPUTE TRANSIENT ONLY'/
» SIMULATION TIME»,19(1H. >, «ENDTIM:',F10.2/
* Tile. STEP' «30(lH.)t 'DT = ',F1C.3/
• MAXIMUM ALLOWABLE TIME STEP •»! 0 < 1H. ) , »DTMAX = ' . F10 .2 /
* MAXIMUM NUMBER OF TIME STEPS' , 12 < 1H .),' HAXNT=' t 11 O/
• TIMC INTEGRATION PARAMETE R' ,1 0 < 1H. ) , 'ALPHA= ' ,F 1 0. 2/
* MAXIMUM ERROR FOR C (JNVERGENCE* , 7 ( 1H . ) , ' tRRMAX= ' »F 10 . 0,1
* MAXIMUM ITERATIONS'»20(1H.> . '1 AX IT= • ,1 1 O/
• TIME STEP CHANGE P ARAMETER • , 1 5 ( 1H. ) , 'CHPARMr • ,F10 .1 )
2001 FORMAT(//lX,72«lh*)//' SOILINtR OUTPUT'//
*lX,72(lH«)//lX,18A1//lX,72(lh*>)
2002 FORMATC' FILE NUMBERS'//
*• GRID INPUT', 10< 1H. ) , • IF GR 0= ', I 1 0 /
*• SOIL PROPERTIES', 8C1H. ) , • IF SOIL= ' »I 1 O/
*' INITIAL CONDITIONS', 7(1H. ), «IFINIT= • ,I10/
*• DUMP FILE'.IOUH.) ,'IFILE = ' ,I10/
*• PSI OUTPUT', 9(1 H. >, MFPOUT=', I10/
*' MOIST OUTPUT', 7 (1H. ), 'IFMOUT='» I10/
*• SOIL PROP OUTPUT', 6UH. ), »IFSOUT=». I10/
*• FLUX OUTPUT', 10(1H. ),'1FLUX = ', I10/
«• ZOJT OUTPUT', 10UH. ), 'IZOUT^', 110)
2001 FORMAT , • NUMNP= « , I lb/
*' NUMBER OF ELEMENTS «COMPUTED> • ,8< 1H . ) , 'NUMEL=« , 11 O/
«• GRADATION PARAM ETER ', 12 < 1 H. ) , • IREGLR-* , 1 1 O/
•' IF IREGLR EQ 1, READ NO. ELEMENTS AND THICKNESS FOR
LAYERS*/
C-14
-------
ooooecoo
00008C10
00008020
00008030
00008040
00008050
00008060
00006C70
00008U80
00008090
00008100
00008110
00008120
00008130
00008143
00008150
00008160
00006170
00008180
00008190
00008200
00008210
00008220
00008230
00008240
00008250
00008260
00008270
00008280
00008290
000083CO
00008310
00008320
00008330
(I n n n ft "*4 n
UUUUOJ^U
00008350
00008360
00008270
fl f) ft fl Q 1 U ft
uu uu o <. o u
00008390
00008400
00008410
00008420
00006430
00008440
00008450
00006460
00008470
00008480
00008490
00008500
00008510
00008520
30008530
00008540
00008550
00008560
00008570
C0008580
00008590
00008600
00008610
*• OTHERylSE, READ NODE. POINT COORDINATES*/
»• WEIGHTED DIFFERENCE OPTION* ,12 C 1M. >, »ISPACE= *» 110 /
*• IF ISPACE EG 1, USc SPECIAL DIf-KLRENCE ALGORITHM*/
»• OTHERylSE, USE STANDARD FINITE DIFFERENCE •>
2006 FORMAT!/* NUMbER OF SPECIAL OUTPUT NODES* , 1 0 ( 1H. ),» NZOUT=
2007 FORMAT!//* SPECIAL OUTPUT NODES*//
*1X,10I10)
2010 FORMAT!/* NUMbER OF SPECIAL OUTPUT TIMES* » 1 01 1H. >,* NOUT=*
2080 FORMAT!/* CONSTANT INITIAL PRESSURE', 4 <1H. ). *PSINIT=* , F10
*, 113)
,110)
.3)
2025 FORMAT!//1X,1011H=>/* SATURATED CONDUCTI V ITY • , 1 0 ( 1H .) , • SATK = * ,
•1PE10.3/1X,20I1H=>//
*• CRITICAL POTENTIAL*, 1011H. ) , *PSICRT = », 1 PE10.3/1 X, 13 1 1H=
2030 FORMAT!//1X,30I1H=>/* MOISTURL RETENTION AND CONDUCTIVITY
•1X,30(1H=)//
«• NJ1&ER OF POINTS ON CUR VE *, : 5 11 H. > , «NCPTS=* , I10//
*• SUCTION MOISTURE RELATIVE CONDUCTIVITY*/
** CL> (-) I-)1/)
2040 FORMAT13F13.3)
2050 FORMAT!//1 DATA FROM INPUT CURVES'//
*« CRITICAL SUCTION PRESSURE*, 13 1 1H. > , «PSI CRT=* ,F 1 0. 3/
*• EFFECTIVE: POROSITY*, 22(in. > ,«poR=*,Fic.4>
2060 FORMAT!//* SPECIAL OUTPUT TIMCS*//
*1X,1P10E12.3>
2069 FORMAT! I10,3F10.3)
2070 FORMAT(2I10,1PE10.3,OPF10.4,OPF10.3J
2071 FORMAT,POKI1) ,DZ!1>
C
URITEUPRT,2000> TIME ,NT, DT,ERR, ITER
URITE(IPRT,2010) ! I ,PSI 1 1 ) ,RMOI STI I ) ,RKi I ),I = 1,NUMNP)
C
IFUFPUUT.LE.O) GOTJ 5
00 4 I=1,NUMNP
PF=O.DO
COP CHANGE NEXT CARD FOR DOUBLE/SINGLE PRECISION
IFIPSKD.LT.O.DO) PF=AL031 0! -PSI ( I) )
WRITEIIFPOUT,2040> I,ZI I) ,PSI 1 I > ,PF
4 CONTINUE
1),
C-15
-------
00008620
00008630
00006610
00008650
00008660
00008670
00008680
00006690
COOOB700
00008710
OC008720
00008730
OOOC&71C
00008750
00008760
00008770
00008780
00003790
00008800
0000881 0
00008820
00008830
00008610
00008850
00008860
00008870
00008880
00008H90
00008900
00008910
00008920
00008930
00008910
00008950
00008960
00008970
00008980
0000 8990
00009000
00009C10
00009020
00009030
00009010
OC00905C
00009060
00009070
00009080
00009090
00009100
00009110
00009120
00009130
00009110
00009150
00009160
00009170
00009180
n n nfl Q 1 Q ft
UUUU 7 A 7 U
00009200
00009210
00009220
00009230
C
5
10
C
20
COP
COP
30
C
35
,RNSTLdN>
IN=IN»1
IF=IP*1
WRITE dFMOUT,2030 > L,Z(IP),R1STL(IN)
CONTINUE
IF drSOUT.LE.O> GOTO 35
I!»=0
IF = 1
DO 30 L=1,NUHEL
IlM=IN+l
RKL03-O.DO
CLOG=O.DO
CHANGE NEXT TWO CARDS FOR DOUBLE/SI NGLL PRECISION
IFIRKL(IN).GT.O.DO) RKLOG =AL051 0 « RKLC IN) )
IF(CHlN) . GT .0.00) CLOG=ALUG1 C(CL( Ii4) >
WRITE ( 1FSOUT.2010 ) 1 ft ,Z ( IP) ,RKL ( IN > ,RKLOS , CL ( IN) , CLOG
IN=IN*1
IP=1P*1
RKL03=O.DO
CLO&-O.DC
CHANGE NEXT TWO CARDS FOR DOUBLE/SINGLE PRECISION
IF(RKLdN).GT.O.DO) RKLOi-ALOGlO(RKLdN))
IF(CLdN).GT.O.DO) CLOG=AL0810 //
*• NODE POTENTIAL MOISTURE K«/)
FORMAT(I6,1X,1PE13.1,1X,OPF12.1,1X,1PE13.1)
FORMATdlO.F 10.3,FiO.D
FORMATdlO»F10.3,F10.7>
FORMAT(I10,OPF10.2,2dPElt;.2,OPF10.3))
END
SUBROUTINE ERRORtA,N,ERRM,ITOL)
C-16
-------
000092*0 C
3000*250 C....
00009260 C
00009270 C
00009Z8Q
OC00929C
000093CO
00009310
00009320
00009130
C0009310 10
00009350
00009360 C
00009370 20
00009380
00009390
00009100 30
00009«tlO
00009120
00009130
00009110 C
00009*60 C
00009170
00009180 C
C00095&C C
00009510 C....
00009520 C
00009530 C
00009510 C
00009550
00009560
00009570 C
00009580 C
00009590 C
00009600
00009610
00009620
00009630
00009610 10
00009650
00009660 C
00009670 C
00009680 C
00009690 20
00009700
00009710
00009720
00009130 30
00009710
00009750
00009760
00009770
00009780
00009790 10
00009800 CBUG
00009810 CBUG
00009P20 CBUG
00009830
00009810
00009850 C
CUMPUTE ERROR IT OL - 0 MAXIMUM ABSOLUTE
ITOL - 1 RMS
DIMENSION A<1>
IFCITOL.EQ. 1 » GOTO 2C
LKRM=O.DO
DO 10 I=1«N
AERR=ABS
SOLVE TRIDIA60NAL MATRIX EQUATION, SEE, FOR EXAMPLE,
PINDER AND GRAY 'FINITE ELEMENT SIMULATION IN
SURFACE AND SUBSURFACE HYDROLOGY*, ACADEMIC PRESS 1977.
DIMENSION A( 3,1), 9( 1) ,C(1) ,W( 1)
IF(IND.EQ.2 ) GOTO 20
REDUCTION REQUIRED IF STIFFNESS IS CHANGING
U(1)=A(2,1)
DO 10 1=2. N
Ihl=I-l
U«I)=A(2,I)-A«1,I)*A(3,IM1)/U(IM1)
CONTINUE
IF(IND.EQ.l) RETURN
SUBSTITUTION
C( l>-b(l)
DO 30 1=2, N
IM1=I-1
C(I)=B(I)-A(l,I)*C(IMl)/y(IMl)
CONTINUE
CiN>=C(N)/W( N)
DO 10 1 = 2, r.
M=N-I»1
MP1=M+1
C«M) = (C
-------
00009870 C
00009880
00009690 C
OOOOS910 C
00009920 C..
00009930 C
00009510
OC00955C
00009960
00009970
00009980
00009990
00010000
00010010
00010020 C
0001Q030
00010010 C
00010C50 C
00010060 C
00010070
00010080
00010090
00010100
00010110
00010120
00010130
00010110
00010150 C
0001016C C
000101 70
00010180
000101SO
0001020D
00010210
00010220 C
00010230 C
00010240
00010250
00010260
00010270
30010280
00010290 C
00010300 C
00010310 C
00010320 C
00010330
00010340
00010350
00010360 CDP
00010370
00010380
00010390
0001 0100
00010110
00010120
00010130
00010140
00010150
00010460 C
00010170 C
SUBROUTINE MASBAHPSI »R MSTL ,OLDML . STARK, DZ. I OUT)
.. COMPJTE AND PRINT MASS BALANCL OF FLOy CALCULATIONS
COMMON/MASS6/STOR,TOTVl,TOTV2,FLUX10tFLUX20
COMMON/ TIMES /TIME.TIME1,NT,ERR.NOUT,TDUT (10) ,^OUT 1, TOUT 1
COMMON/DEVICE/ IRQ «IPRTt IFILE
COM10N/INFO/NUMN?,NUMEL.NPi' ,OLDML< 1) ,DZ( 1 )
IF(IOUT.NE.-l) GOTO 5
INITIALIZE TERMS
TuTVl^O.DO
TOTV2=O.DO
STOR=O.DO
FLUX 10 = STARK«1)«(1.DO + (PSI(2>-P SKI))/ DZ<1»
IF(OZd).oT.O.DO) FLUX10 = -FLUX10
FLUX20 = STARK (NUKEL)» ( 1 . D0*( PS I < NUM^P) -PS I < NUMtlL) >/OZ( MUCE.L) )
IF(DZ (NUHEL) .LT.0.00) FLUX2 0 = -FLUX 20
RLTU3N
TOP FLUX
5 FLUX1=3TARK(1)*(1.DO+(PSK2)-PSI(1))/DZ(1))
IF(DZ (D.ST.O.DO) FLUXl=-FLUxl
V OL1 =( FLUX 1 • ALPHA +FLUX10* AM 1)*DT
FLUX1'"> = FLUX1
TOTV1=TOTV1»VOL1
BOTTOM FLUX
FLUX2=STARK« NUMEL)*( 1 .0 0+ CPSHMUMNP )-PSI ( NUM£.L> > /DZ -OLDML ( INI ) +RMSTL t IN2)-OLDML (IN2)) *DEL2/
* 2. DC
10 CONTINUL
STOR=STOR»DSTOR
IFtIOUT.NE.1) RETURN
RATEST=DSTOR/DT
COMPUTE FLUX AND VOLUME ERRORS
C-18
-------
00010480 C
00010490 COP
0001050C
00010510
C0010520
00010530
00010540 C
0001055C
CHANGE NEXT THREE CARDS FOR
EFLUX=ABS(-FLUX1«FLUX2+RATEST
£VOL=AbS(OSTOR-VOLl+VOL2>
ETuT=AbS«STOR-TOTVl*TCTV2 )
E»EL=ETQT/STOK
URITEUPRT.2020) FLUX1.FLUX2.
00010560 * OSTOR.EVOL.T
00010570 C
00010580
00010590 C
COOlOfcCf 2020
00010610
00010620
00010630
00010640
00010650
00010660
00010670
00010680
00010690
00010700
00010710
00010720
00010730
00010740
00010750
00010760
30010770
00010780
00010790 C
DCOlOfalC
00010830 C
00010640 C....
00010650 C
00010660
00010870
0001088C «
00010690 •
00010900
00010910 19
00010920
00010S30 C
00010540
00010S50
0001 0960
00010970 1020
00010S80 2020
RETURN
DOUBLE/SINGLE PRECISION
)
FATEST.LFLUX.VOL1 .VOL2 t
OTVl,TOTV2,STORȣTl;T,EREL
FORMAT«//» VOLUME BALANCE CALCULATIONS*//
« TOP FLUX «,!PEli.4/
• BOTTOM FLUX».lPE12.t/
• STORAGE RATE*,lFri2.W
20X.12I1H-)/
• tRROk '.1PE12.4/
• VOLUME IM '.1PE12.4/
• VOLUME OJT «,lPt_12.4/
• STORAGE VOLUME «.1PE12.4/
20X»12<1H-)/
//
• ERROR «.1PE12.4///
* CUMULATIVE CHANGES'//
* VOLUME IN (-) «,1PL12.4/
' VOLUME OUT »tl?E12.4/
• STORAGE «flPcl2.4/
20X.12C 1H-)/
' ERROR »,1PE12.4//
' RELATIVE ERROR '.OPF12.6)
END
SUBROUTINE NEUBC(PSI)
READ AND CHANGE PRESSURE BOUNDARY CONDITIONS
DIMENSION PSK1)
COMMON/INFC/NUMN'.NUMEL.NPMl,
NPH2.NUMEL2,AM1,ISPACE ,ITOL,
DEPTH .ENDTIM.DT «D TMAX » ALPHA.P SIBOT .H » NZ OUT.
PSINIT.NCPTS,cRRMAX,MhXIT.LHPARM,MAXNT.ISTLDY
COMHOM/DEVICE/IRD.IPRT, IFILE
REAO< IKD,1??0) H.PSIBOT
klrtlTE (IPHT,^020) H.PSIfiOT
APPLY PRESSURE BOUNDARY COND
PSH1 ) = H
PSI(NUMNP)=PSIBOT
RETURN
FORMAT(3F10. C)
ITIONS
FOKMATl//lX,30llHr)/t BOUNDARY AND INITIAL COND 1 T 10 NS • /IX .3 0 ( IHr ) /
C0010590 */• HEAD IN I MPOUNDMENT* ,2 0< 1H
00011000 «
00011C10
C0011020 C
00011 040
00011060 C....
00011070
00011080
00011C90
.), «H='.F10.2/
• UNDERLYING SOIL SUCTION PRE t>SURE • .5 ( 1H. ) , 'PSIB OT= •» Fl 0. 3>
END
SUBROUTINE V SCALE C A.D »'J ,C>
MULTIPLY VECTOR A 3Y SCALAR B
DIMENSION All)tC(l)
IHN.LE.O) RETURN
DO 10 1=1, N
AND RETURN IU C
C-19
-------
00011100
00011110
00011120
00011130
0001 1150
00011170
0001118C
0001 1190
0001 1200
00011210
00011220
C001 1230
000112*0
00011260
G001 1280
00011290
0001 1300
0001131Q
00011320
00011330
000113*0
00011350
00011370
0001 139C
00011 tOO
00011*10
00011*20
00011*33
0001 1**0
00011*50
00011460
0001 1*80
0001 1500
00011510
00011520
00011530
00011510
00011550
00011560
00011570
00011590
00011610
00011620
00011630
000116*0
0001 1650
00011660
00011670
0001 1680
0001 Ifc90
00011700
00011710
C( I)=A( I>*6
10 CONTINUE
RETURN
ENO
SUBROUTINE CLEAR(A.N)
C.... FILL A WITH ZEROS
DIMENSION A(l)
IK(N.LE.O) RETURN
DJ 13 1=1, N
M I)=0.00
10 CONTINUE
RLTURN
END
SUBROUTINE HADD
C.... PROGRAM TO RETURN THt SUM OF A AND B IK C
DIMENSION A< 1),B< 1) ,C<1 >
IF(N.LE.O) RETURN
DO 10 I-1»N
C( D-AC I) + B« I)
10 CONTINUE
RCTUHN
ENU
SUBROUTINE COFYtA,N,B)
C.... COPY VECTOR A INTO B
DIMENSION A(l), 6(1)
IF(N.LE.O) RETURN
DU 10 1=1, N
B(I) = Ad )
10 CONTINUE
RLTURN
END
SUBROUTINE DUMP1 ( NAME ,A RR AY *N )
C.... PRINT ARRAY
COMMON/DEVICE/IRD,IPRT,IFILE
DIMENSION ARRAYtl ),NAME«2)
URITEI IFILE, 20CO) NAME
IF(N.LE.O) KETURN
UKITECIFILE,201D) ( ARRAY( I) , I =1 , N)
RETURN
2000 FORMAT«//« DUMPING ARRAY «,2A4,« . . .»>
2010 FORMAT(lP8E:i0.31
END
C
C-20
-------
00011720
00011730
00011740
00011750
00011760
00011770
00011780
00011790
00011600
OOOlldlC
00011620
00011830
0001181C
00011850
00011860
00011670
0001 1880
0001 1890
00011900
00011910
SUBROUTINE ZOUTtlZOUT.PSI.RMOIST)
C
C.... URITE PSI AND RMOIST AT SPCCI6L NODES LACh TIML STEP
C
COMMON/FILES/IFGRO,IFSCIL,IFINIT»IFPOUT,IFMOUT,IFSOUT, IFLUXf IFZOUT
COMMON/INFO/NUMNP»NUMt.L.NrMl»NPM2,NUKEL2,AMl,ISPAC£,ITOL»
* DEPTH. ENDTIM.DT.DTM AX » ALPHA«PSIBOT »H »f;Z OUT,
* PSINIT,NCPTS,ERRHAX»MAXIT,CHPARM,MAXNT, ISTEDY
COMMON/TIMES/TIHE.TIMEl,NT,ERR.NOUT.TOUr(10),NOUTl,TOUTl
COMMON/ERRORS/IERRtlHARN
D I MEM SI ON IZOUTtl ),PSH1),RMOIST(1)
00 10 I=1,UZOUT
N=IZOUT(I>
URITE( IFZOUT .20 00) N, TIME»PSI (M ) tRMOI ST ( N )
10 CONTINUE
RE.TURN
2000 FORMA T( 110, 1 PE10. 3, OPF1 0. 2» OPFlO.t)
END
C-21
-------
APPENDIX D
EXAMPLE INPUT AND OUTPUT FOR SOILINER MODEL
D-l
-------
09/30/83 10:04:21 3CAGD.TEST3A.OUTPUT.DATA
•A**********************************************************************
SOILINER OUTPUT
ft***********************************************************************
TEST3A 50 1 CM ELEMENTS, 51 NUMNP YOLO INFILTRATION H=25.
•ft**********************************************************************
FILE NUMBERS
&KID INPUT ...... ,...IFGRD= 21
SOIL PROPERTIES ........ IFSOIL= 22
INITIAL CONDITIONS ....... IFINIT= 23
DUMP FILE .......... IFILE= 10
PSI OUTPUT.... ..... IFPOUT= 31
HOIST OUTPUT ....... IFMOUT= 32
SOIL PROP OUTPUT ...... IFSOUT= 33
FLUX OUTPUT .......... IFLJX= 3T OUTPUT .......... IZOJT= 0
TEMPORAL DISCRETIZATION PARAMETERS
INITIAL CONDITION 3A^M.... ...... INITAL= D
IF INITAL EQ 1, ^EAD INITIAL PSI FOR EACH NODE
OTHERWISE, READ CONSTANT INITIAL PSI
STEADY STATE PARM ............ ISTEDY= 0
IF ISTEDY EQ 1, COMPUTE STEADY STATE SOLUTION
OTHERWISE:, COMPJTE TRAMSIENT OMLY
SIMULATION TIME ................... ENDTIM = 200000.00
TIME STEP .............................. DT= 0.100
MAXIMUM ALLOWABLE TIME STEP .......... DTMAX= 1000.00
MAXIMUM NUMBER aF TIME STEPS ............ MAXNTs 1000
IIME INTEGRATION PARAMETER .......... AL3HA= 0.50
MAXIMUM tRROR FUR CONVERGENCE ....... ERRMAXr 0.1000
MAXIMUM ITERATIONS .................... MAXIT= 10
TIME STEP CHANGE PARAMETER ............... CHPARM= 25.0000
NUMBER OF SPECIAL OUTPJT TIMES .......... »JOJT= 5
SPECIAL OUTPUT TIMES
1.000F.+ 03 l.OOOE+Of t.OOOE+0* 1.000E»05 2.000E*Ob
SPATIAL DISCRETIZATION PARAMATEKS
NUMBER OF NODE POINTS NUMNf= 51
NUMBER OF ELEMENTS (COMPUTED) NUMEL= 50
6RADATION PARAMETER IREGL1'= 1
IF IREGL* EQ 1, ^EAD >JD. ELEMENTS AND THICKNESS FOR .AYERS
OTHERylSE, READ NODE POINT COORDINATES
WEIGHTED DIFFERENCE OPTION ISPACE= 0
IF ISPACE EQ 1, JS£ SPECIAL DIFFEhENCE ALGORITHM
OTHERWISE, USE STANDARD FINITE DIFFERENCE
D-2
-------
NUMBER OF SPECIAL OUTPUT NODES WGUT =
CRIO DA1A
ELMENT
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
NODE Z
1 0.0
-1.000
2 -1.000
-1. 000
3 -2.000
-1.000
4 -3.000
-1.000
t> -4.000
-1.000
6 -5.003
-1.000
7 -6.000
-1.000
8 -7.000
-1.000
9 -8.000
-1.000
10 -9.000
-1.000
11 -10.000
-l.OUO
12 -11.000
-1.000
13 -12.000
-1.000
14 -13.000
-1.000
15 -14.000
-1.000
Ifc -15.000
-1.000
17 -16.000
-l.OUO
It* -17.000
-1.000
19 -18.000
-1.000
20 -19.000
-l.OUO
21 -20.000
-1.000
22 -21.000
-1.000
23 -22.000
-1. 000
24 -23.003
-1.000
25 -24.000
-1.000
26 -25.000
-1. 000
27 -26. ODD
-l.OUO
28 -27.300
DZ
-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
DZN
1 .000
1 .DOO
l.COC
1.000
1 .300
1.000
1 .000
l.DOO
1 .000
1.000
1 .OOP
l.OCO
1.000
1 .000
1.000
1.000
1.0 CO
1.000
1 .000
1 .000
1.000
1 .000
1 .0 00
1.000
1.000
1.000
1.000
D-3
-------
28
29
30
31
32
33
31
35
36
37
38
39
4 D
41
42
43
44
45
46
4 7
48
49
50
SOIL
-1.000
29 -28.000
-1.000
3U -29.000
-1.000
31 -30. COD
-l.OCO
32 -31.003
-1. 000
33 -32. COD
-1 .000
34 -33.000
-1.000
35 -34. ODD
-1. 000
3fc -35.033
-1.000
37 -36.000
-1. 000
38 -37.000
-1.000
39 -38.003
-l.COO
40 -39.003
-1.000
41 -40.000
-1.000
4^ -41.C03
-1.000
43 -42.000
-1.000
44 -43.003
-l.COO
4b -44. ODD
-l.OCO
46 -45.000
-1.000
47 -46.000
-1.000
48 -47. COO
-1 . 000
49 -48.000
-l.COO
50 -49. COO
-1. 000
51 -5C.COD
PROPERTIES
ELEMENT SATK
1 1
2 1
3 1
4 1
5 1
fc 1
7 1
8 1
1 I
10 1
-1.000
-1 . OOC
-1. 000
-1 . 000
-1.000
-1. 000
-1 . 000
-i . oca
-1 . OOC
-1 . 000
-1.000
-1.000
-1 . 000
M. 000
-1 . 000
-1.000
-1. 000
-l.OOC
-1 . 000
-1. 000
-1. 000
-1 . 000
FOR PS
1.23DE-05
1 .230L-05
1.230E-05
1.230o-05
1.230L-05
1 .233E.-C5
1 .233E-05
1 .23Ct>05
1.250L-05
1.230c>05
1.000
1.100
1 .300
l.OQO
1.000
1 .000
1.000
1.000
1 .003
1 .000
1 .000
1.000
l.OOC
l.OCO
1 .000
1.000
1.000
1 .000
1 .0 PO
l.COO
1 .0 CO
1 .300
ICRT
0.4950
0.4950
0.4950
0.4950
C.4950
0.4950
0.495C
0.4950
0.495C
0.4950
-1.000
-l.OOC
-1 .000
-1.000
-l.COO
-1. OOC
-1.000
-1.000
-l.'OOD
-l.'OOO
D-4
-------
11
12
13
16
17
18
19
20
21
22
23
21
25
26
27
28
2^
3C
31
32
33
34
35
36
37
38
3V
40
41
42
43
44
45
46
47
4P
49
50
1 1.239E-05
1 1.233E-05
1.230L-05
1.230E-05
1.230E-05
1 . 2 3 3 E -0 5
1.233E-05
1.233E-05
1.230E-05
1.230L-05
1 1.230L-05
1 1.230E-05
1.233E-05
1.233E-05
1.230L-05
1.230E-05
1 1.230C-05
1 1.230E-05
1 1.230£-05
1 1.230E-05
1 1.230E-05
1 1.230E-05
1 1.230E-05
1.23DE-05
1.23DE-05
1.233E-05
1 .230E-05
1 1.23CE-G5
1 1.230c:-05
1 1.230E-05
1.230E-05
1.230E-05
1.230E-05
1.230t-05
1.230E-05
1.230E-05
1 1.230E-05
1 1.230E-05
1 1.230E-05
1 1.230E-05
0.4950
0.4950
0.495D
0.4950
0.4950
0.4950
0.4550
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.495C
C.4950
0.4350
0.4950
0.4950
0.4950
0.4950
0.4950
3.4950
0.4950
0.4950
0.4950
0.4950
0.4P50
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
-1.000
-1.000
-l.OOC
-1.000
-1.000
-l.OOC
-1.000
-l.OOC
-1.000
-1.000
-1.000
-l.COO
-1.000
-1.000
-1.000
-l.OOC
-1.000
-l.OOC
-1.OOQ
-l.OOC
-1.000
-WOOD
-1.000
-l.OOC
-1.000
-1.000
-1.000
-1.000
-l.'OOO
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-1.000
-l.OOC
-1.000
-1.000
CONSTANT INITIAL P^ESSURE... . PSINIT- -600.000
BOUNDARY AND INITIAL CONDITIONS
HEAD IN IMPOUNDMENT H- 25.00
UNDER-LYING SOIL SUCTION PRESSURE PSISOT= -600.000
TIME= 0.0
ITER= 0
TIME STEP=
0 DT= l.OOOE-01 ERR=
0.0
NODL
POTENTIAL
MOISTURE
D-5
-------
1
1
3
4
5
6
7
8
9
1C
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
3C
31
32
33
31
35
36
37
38
39
40
41
42
43
44
45
46
47
48
.49
50
51
**»*.**
2.5300E+31
-6.00COE+02
-6.0GOOE+02
-6.0000E+C2
-6.0COOE+C2
-6.0030E+D2
-6.0000E+02
-6.0000L+02
-6.0000E+C2
-6.00001+02
-6.DCCOE+02
-6.0000E+02
-6.0UOOE+02
-6.0COOE+02
-6.0000E+92
-6.QCOOE+32
-6.00001+02
-6.0000E+02
-6.0000E+C2
-6.0000E>D2
-6.0003E+32
-6.00COE+02
-6.0000E+C2
-6.0000E+02
-6.00001+02
-6.0000E+02
-6.0000E+02
-6.0000E+02
-6.0000E+02
-6.00001+32
-6.0000E+C2
-6.00001+02
-6.0000E+02
-6.0000E+02
-6.3000E+02
-6.0000E+02
-6.0000£+n2
-6.00CCE+02
-6.0000E+02
-fa.OOOOE+02
-6.0000E+C2
-t.OOOOE+02
-6.0000E+02
-6.0000E+02
-&.OOOOE+02
-6.00CO£+02
-6.00COE+02
-6.000QE+02
-6.0000E+02
-6.0000E+02
-6.0000E+02
********************
D.4950
0.2376
0.2376
0.2376
0.2376
3.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.237b
3.2376
3.2376
0.23/6
0.2376
0.23/6
0.2376
0.2376
0.2376
0.2376
0.2376
3.2376
0.2376
0.2376
0.2376
P.2376
0.2376
0.2376
0.2376
0.2376
0.2376
3.2376
0.2376
C.2376
0.2376
3.2376
0.2376
0.2376
0.2376
0.2376
3.2376
3.2376
0.2376
0.2376
0.2376
3.2376
C.2376
3.2376
*************
1 .2300E-35
1.8511E-08
1.8511L-08
1.8511E-08
1.8511E-08
1.S511E-03
1.8511E-08
1.8511E-08
1.8511E-08
1.B511E-01*
1.8511E-03
1.8511E-OB
1.8511E-08
1.8511E-06
1.S511E-D8
1.8511E-08
1.8511E-08
1.8511£-08
1.8511E-06
1.8511E-OS
1.8511E-08
1.8511E-08
1.8511E-OU
1.8511E-08
1.3511E-33
1 .8511E-08
1.8511E-OS
1.8511E-08
1.8511E-08
1.8511E-33
1.8511E-06
1.8511E-08
1.8511E-OS
1.8511E-OB
1.8511L-09
1.8511E-08
1.8511E-08
1.8511E-08
1.3511E-38
1.8511E-39
1.8511C-OH
1 .6511L-08
1.8511E-08
1 .8511E-03
1.8511E-08
1.8511E-08
1.8511E-03
1.8511E-08
1.8511E-03
1.8511E-08
1.8511E-OB
********************************
TIME= 1.000E+C3 TIME STEP=
\
ITER= 2
35 DT=
3.4b2E+01
1.724E-02
NODE
POTENTIAL
HOISTURE
D-6
-------
1
2
3
i»
5
6
7
b
1C
11
12
1 3
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
3 7
38
39
40
41
42
43
44
4i>
46
47
48
49
50
51
2.5000-: + Cl
-1.1466E+P1
-2.13i2E+C2
-5.5319E+02
-5.9774£ + C 2
-5.9991E+02
-fc.OPOOS+02
-6.00CO£+C2
-6.0000E+02
-f.OOOOE+02
-6.3COOE+02
-6.0CCOL+02
-6.0000E+02
-6.0 CODE + 02
-6.0000E+02
-6.00COE+02
-6.0COOE+02
-6. OOOOE>02
-6.0000E+02
-6.0000i+02
-6.0000^+02
-6. 0000:1 + 02
-6.0000E+02
-f.OOOOE+02
-6.0000:1 + 02
-6.0000i+02
-6.0000-1+02
-6.0000£+02
-6.000CE+02
-6.0COOE+02
-t .OOOOE + 02
-£.OOOOE: + Ci2
-6. 00 00 £ + 02
-6.00COE+P2
-6.0000E+C2
-6.0000E+C2
-fa.OOOOE: + 02
-6.0000E+C2
-6.0 G CO E* 02
-6.0CCO£+02
-6.00CO£+02
-6.0000E+02
-6.0000E+02
-6.0000L+C2
-6.00COE+C2
-f, *OOCO£ + 02
-6.0CCCE+C2
-6.00COE+C2
-6.0000£+02
0 .49=30
0.4780
0.2991
0.2417
C.2378
0.2376
0.2376
0.2376
0.2376
0.2376
C.2376
0.2376
D.2376
D.2376
0.2376
0.2376
0.2376
0.2376
0.2376
C.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
C.2376
C.237b
0.2376
0.2376
0.2376
0.23.76
0.2376
D.2376
0.2376
0.2376
C.2376
0.2376
D.2376
0.2376
C.2376
0.2376
C .2376
0.237b
0.2376
0.2376
0.2376
0.2376
1.2300E-05
7.6775E-OS
1.1473L-07
2.1 369£-0:)
1.8635E-06
1.8516E-&8
1.6512E-D8
1.8512E-08
1 .8511£-0 3
1 .8511E-03
1.85UE-08
1.65 HE -Oa
1.8511L-38
1.8511E-08
1.8511E-OS
1.8511E-Oa
1.8511E-03
1.&511E-OS
1.8511E-08
1.8511E-08
1.8511E-08
1.8511E-08
1.35HE-03
1.8511E-OS
1.8511E-OS
1.8511E-08
1.8511E-06
1.8511E-08
1 .8511E-06
1.8511E-08
1.8511E-08
1.8511E-08
1.6511L-08
1.8511E-08
1.8511E-08
1.8511E-08
1.8511E-08
1.8511E-09
1.8511E-OS
1.8511t-OS
1. 8511006
1.8511E-03
1.8511E-08
1 .8511E-08
1.8511E-08
1.8511E-08
1.8511E-05
1.8511E-08
1.8511E-08
1.8511E-08
1.8511E-08
VOLUME BALANCE CALCULATIONS
TOP FLUX 3.6408E-G4
BOTTOM FLUX -l.fi511E-08
STORAGE RATE 3.62S1L-0*
ERROR
D-7
-------
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
1.2547E-02
-i.406SE-C7
1.2553E-02
5.0775L-06
CUMULATIVE CHANGES
VOLUME IN
VOLUME OU
tTORAGF
ERROR
RELATIVE
TIME = 1
1TER =
NODE
1
£.
3
4
5
6
7
8
q
10
11
12
13
It
15
16
17
18
19
20
21
22
23
24
25
26
27
26
29
30
31
32
33
(-) 3.0680E-01
T -1.8511E-05
3.0623E-01
5.8619E-04
ERROR 0.301911
.OOOE+04 TIME STEP=
2
POTENTIAL HOI
2.5000:1 + 01
2.0385E+01
1.1389E+01
4.0702E+00
-3.1316E+00
-1.2718E+31
-3.2157E+01
-9.2790E+01
-2.7631£+C2
-4.949JE+02
-5.791£>E: + 02
-5.9658E+02
-5.9948E+G2
-5.9992E+02
-5.9553E+02
-6.0COOi+C2
-6.00COE+02
-6.0000E+02
-6.0000E+02
-6.0000E+02
-b.OOOOE+02
-6.00001+02
-6.0000E+02
-6.0000E+02
-b. 0000-1 + 02
-fa.OLCOE+02
-6.0COOE+C2
-6.0000E+C2
-6.0000E+02
-6.0000i+02
-6.0COCE+02
-6.C3COE+02
-6.0000E+02
109 DT=
STURE
D.t950
O.«950
C.4950
0.4950
0.494 1
0.4751
0.4341
3.3603
3.2816
3.2475
3.2394
C.2379
C.2376
0.2376
3.2376
0.2376
C.2376
0.2376
0.2376
3.2376
3.2376
0.2376
0.2376
0.2376
3.2376
3.2376
3.2376
0.2376
D. 2376
3.2376
0.2376
C .2 376
0.2376
1.284E+02 E^R= 1.969E-02
K
1.2300E-D5
1.2300E-05
1.2300E-05
1.2300E-05
1.1598E-05
7.1375E-06
2.5973E-06.
4.8471E-07
7.2709E-OS
2.6011E-08
1 .9705E-03
1.8699E-08
1.8540E-OB
1 .6516E-08
1.8512E-OS
1.8512L-33
1.8512E-08
1.8512E-08
1.6512E-08
1.35X1E-08
1.8511E-08
1.8511E-08
1.8511E-38
1.8511E-Oi
1.8511E-03
1.8511E-08
1 .8511E-08
1.8511E-06
1.8511E-09
1.3511E-D3
1 .8511E-3S
1.8511E-09
1.8511E-08
D-8
-------
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
-6.0000E+02
-6.0000i+02
-6.0000i+02
-6.0000i+02
-6.0000E+02
-6.0000E+02
-6.0000i+02
-6.0000E+02
-6.0000E+02
-6.0000E+G2
-6.00001+02
-6.00001+02
-6.00COE+02
-6.0000E+C2
-e.oocoi+02
-6.0000i+02
-6*OOOOE+02
-6.00001+02
3.2376
0.2376
0.2376
C.2376
C.2376
0.2376
0.2376
0.2376
0.2376
0.2376
3.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
1.3511E-08
1.8511E-08
1.8511E-08
1.85UL-08
1.8511E-08
1.8511E-03
1.6511E-09
1.8511E-OB
1.8511E-Oa
1.8511E-08
1.3511L-09
1 .8511E-03
1.8511L-OB
1.8511E-08
1.85ME-Q8
l.a511E-03
1.8511E-08
1.8511E-08
VOLUME BALANCE CALCULATIONS
IOP FLUX 6.9067E-05
BOTTOM FLUX -1.8511E-08
STORAGE RATE 9.8334E-05
ERROR
2.9248E-05
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
1.2635E-02
-2.3771E-06
1.2627E-02
9.62S5E-06
CUMULATIVE CHANGES
VOLUME IN <-)
WOLUME OUT
STORAGE
ERROR
RELATIVE ERROR
1.6463E+00
-1.8512E-04
1.641SE+00
4.8181E-03
0.002935
TIME= 4.000E+04 TIME STEP= 175 DT= 5.115E+02
1TER= 2
******************»****************»tlkA>AAtAtttAI
NODE
1
2
3
4
5
POTENTIAL
2.50COi+01
2.4162E+01
1.6944E+01
1.54021*01
1.2220E+01
MOISTURE
0.4950
0.4950
3.4950
3.4950
0.4950
K
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-35
1.2300E-05
3.534E-02
>********
D-9
-------
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
3b
37
38
35
4 3
41
42
43
44
45
tfe
47
48
49
50
51
9.0307E+00
5.S390E+CO
2.6502-1 + 00
-5.5438E-01
-3.9420E+00
-7.9113E+00
-1.3035E+D1
-2.0385E+01
-3.2131L+01
-5.3065E+01
-9.3773E+C1
-1.7276E+02
-2. 99321 + 02
-4.3570E+02
-5.2863E+32
-5.7326E+02
-5.9078E+02
-5.9698E+02
-5.9904E+C2
-5.9970E+02
-5.9991E+02
-5.9997E+P2
-5.9999E+02
-5.9999i+02
-5.9999E+02
-5.9999E+02
-5.999iE+02
-6.000DE+02
-6.0000E+02
-6.03COE+C2
-6.0000E+02
-6.0000E+02
-6.0000:+02
-6.0000E+02
-6.COOCE+C2
-6.000QE+02
-6.0000E+02
-6.0COOE+02
-6.0000E+02
-fe.OOCOE+02
-b.OOOCE+02
-6.0000E+D2
-6.00CQE+02
-6.0000E+02
-6.0000E+02
-6.0000E+C2
0.4950
0.4950
0.4^50
P. 4 950
C.4932
0.4860
3. 4744
3.4577
0.4342
0.4016
0.3595
0.3139
0.2767
0.2544
0 .2440
0.2399
3.2384
0.2378
3.2377
C.2376
0.2376
3.2376
0.2376
0.2376
0.2376
0.2376
0 . 2376
C.2376
0.2376
0.2376
0.2376
3.2376
0.2376
0.2376
0 . 237fe
C.2376
3.2376
C.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
0.2376
3.2376
1.2300E-05
1.2300E-05
1.2300E-OD
1 .2300E-05
1.1274E-05
9.3738E-06
7.0066E-06
4.6117E-0&
2.6003E-OS
1.2220E-05,
4.7610E-07
1.6568E-07
6.3158E-03
3.2577E-03
2.3154E-08
2.0065E-06
1.9C25E-08
1 .BS73E-03
1. 8564[>08
1.8528E-08
1.8517E-08
1.8513E-08
1.6512C-03
1.8512E-OS
1.8512E-08
1.6512E-03
1.8512E-03
1.8512E-G8
1.8512E-08
1.8512E-06
1.8512E-08
1 .8512E-03
1.8512E-D8
1.&512E-OB
1.6512£-08
1.8512E-08
1.8512E-D3
1.8512E-08
1.8512E-03
1.8512E-08
1.8512E-08
1 .S512L-03
1.8512E-08
1.8512E-08
l.&511E-Oa
1.8511E-08
VOLUME BALANCE CALCULATIONS
TOP FLUX 2.26C6E-05
BOTTOM FLUX -1.8516L-08
STORAGE RATE 5.1723E-05
ERROR
2.9095E-05
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
2.6478E-02
-9.4717E-06
2.6457E-02
3.0821E-05
D-10
-------
CUMULATIVE CHANGES
VOLUME IN (-)
VOLUKE OUT
STORAGE
ERROR
RELATIVE ERROR
3.6517E»00
-7.4G53t>04
3.644SK>00
7.8416E-U3
0.002152
r*********************************************
TIME= 1.000E+C5 TIME STEP= 238 DT= 2.031E+01 ERR = 7.258E-02
ITER= 2
***********************************************************************
NODE
3
4
5
6
7
8
9
10
11
1?
13
It
15
16
17
18
19
2C
21
22
23
2*
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
POTENTIAL
2.5000:1*01
2.07121+01
2.08P4E+01
1.9219E*01
1.7275E+01
1.5338£+01
1.34031*01
1.1466E+01
9.5441E+00
7.6149E+00
5.6794£+00
3.7810E+CO
2.0092E+00
3.7577E-C1
-1.24741 + 00
-3.0354i+00
-5.2093E+00
-7.6458E+GO
-1.0525E+01
-1»4082E»01
-1.8&67:+31
-2.4849E + 01
-3.3585E+&1
-4.&520E+01
-6.6485E+01
-9.8079E+01
-1.47571+02
-2.1994E+02
-3.1167E+02
-4.06361+32
-4.8436i+02
-5.3710E+02
-5.6795E+02
-5.8437E+P2
-5.9260E+02
-5.96581+02
-5.9844i+02
-5.9930E+02
-5.9969E+C2
MOISTURE
D.4950
0.4950
0 . 4 95 0
C. 4950
0.4950
0.4950
0.4950
D . 4 95 0
0.4950
0.4950
0.4950
0.4950
0 . 4 95 0
0.4950
0.4950
0.4942
0.4913
0.4U02
3.4720
0.4615
0.4482
0.4315
0.4107
3.3853
C.3561
0.3254
0.2970
0.2741
0.2533
0.2486
0.2432
0.2*03
0.2389
0.23S2
0.2379
D.2377
j.2377
C.2376
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.23COE-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2156E-05
1.1634E-35
1.0705E-05
9..50b4£-0b
8.1061E-06
6.5910E-06
5.0689L-05
3.6531E-D&
2.4433E-DS
1.5034E-06
•>.4761£-07
4.4104E-07
2.1804E-07
1.0S57E-07
5.8&15E-08
3.6842E-03
2.7C22E-08
2.25.13E-03
2.0397E-08
1.9395E-08
1.8922E-03
1.6 /OOE-08
1.8597E-08
1.6550E-08
1.8529E-08
D-ll
-------
40
4 1
42
43
44
45
4£
47
48
49
50
51
-5.99661+02
-5.9994E+C2
-5.9997E+02
-5.9998E+02
-5.9599E+02
-5.99991+32
-5.9999E+C2
-6.0000E+02
-6.0000E+02
-6.00001+02
-fa.0003£+02
-6.0000E+02
0.2376
C.2376
0.2376
0.2376
D.2376
C.2376
0.2376
0.2376
0.2376
0.2376
3.2376
0.2376
1.8519E-03
1.8515E-CS
1.6513c>08
1.8512E-08
1.8512E-08
1.8512E-33
1.8512E-OS
1..6512E-D8
1.0512E-06
1.3512E-08
1 .8512E-08
1.B511E-08
VOLUME BALANCE CftLCULATIONS
TOP FLUX 6.5047E-05
BOTTOM FLUX -1.8534E-08
STORAGE RATE 3.5356E-05
ERROR
2.9709E-05
¥OLUKE IN
VOLUME OUT
STORAGE VOLUME
ERROR
7.1503E-04
-3.7647E-07
7.18181-04
2.7665E-06
CUMULATIVE CHAN3ES
VOLUME IN t->
VOLUME OUT
STORAGE
ERROR
RELATIVE ERROR
6.1434E+00
-1.8519E-03
6.1338E+OC
1.1426E-02
0.001863
****«***«»•*************»*****»****»*****»****»*»*****«
TIME= 2.000E+05 TIME STEP= 340 DT= 3.5&7E+02
ITER= 2
*******••»**»*«
9.142E-02
NODE
1
O
C.
3
4
5
6
7
8
V
10
11
POTENTIAL
2.5000E+01
2.1246E+01
2.1872E+01
2.0621E+01
1.9411E+01
1.8008E+C1
1.6607E+01
1.5204E+01
1.3&1&1+0 1
1.2421i+01
1.1019E+01
MOISTJRE
0.4950
0 . 4 95 0
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
0.4950
3.4950
K
1.2300E-05
1.2300L-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1 .2300E-05
1.2300E-05
1.2300t-05
D-12
-------
12
13
11
Ib
16
17
18
19
20
21
22
23
21
25
26
27
28
29
30
31
32
33
31
35
36
37
3b
39
10
11
12
13
11
15
16
17
18
19
50
51
9.6515E+ 00
8.11611+CO
7.31fa8i+30
fa.2892i+00
5.3215E+00
1.2836£+00
3.2181E+00
2.2666i+00
1.2812E+00
2.6tt01i-01
-7.3933E-01
-1.75681+00
-2.9237£+00
-1.2899E+00
-5.8026E+CC
-7.1711E+00
-9.3&8&1+00
-1.1567i+31
-1.11731+01
-1 .7336E+01
-2.1269E+01
-2.62851+31
-3.2851E+01
-1.16721+01
-5.3815E+01
-7.0856E+01
-9.5311-: + Cl
-1.2935i+C2
-1 .7557E+02
-2.3538E+02
-3.0535E+02
-3.77931+02
-1*1393i+02
-1.9691i+02
-5.3531E+G2
-5.6109E+02
-5.7751i+:2
-5.8790i+02
-5.9171i+P2
-6.0000E+02
0.1950
?. 195C
D.I 950
0.1950
3 . 1 95 0
0.1950
0.1 950
0.1950
0.1950
3.1950
0.1950
3.1919
0.1913
0.1928
0.1903
0.1b70
3.1828
3.1778
0.1717
3 .1615
0.1558
0.1153
3.1329
0.1180
3.1006
0.3836
3.3535
3.3351
3.3128
0.2921
0.2751
3.2625
3.2533
0.2172
0.213-3
0.2109
0.2335
0.2386
3.2330
0.2376
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300E-05
1.2300£-0:>
1 .2300E-05
1.2300E-G5
1.2300E-05
1.2300E-05
1.2038E-05
1.1671 £-05
1.1125E-05
1.0121E-05
9..5923E-06
8.6555£>0b
7.6329E-0=>
b.5559E-Oi
5.1630E-3S
1.3972E-0&
3.1031E-06
2.5208E-06
1.7800E-06
1.1919t-0:>
7.6286E-07
1.D556E-07
2.7515E-07
1.6107E-07
9..6377E-08
6.0977E-OB
1.1673E-03
3.1519E-03
2.5826E-08
2.2611E-08
2.0810E-08
1 .9b03E-Oi
1 .9190E-08
1.8802E-08
1.8511E-08
VOLUME BALANCE CALCULATIONS
TOP FLUX 5.8178E-05
BOTTOM FLUX -1.1682E-07
STORAGE RATE 2.7076E-05
ERROR
3.1513E-05
VOLUME IN
VOLUME OUT
STORAGE VOLUME
ERROR
9.69S3E-03
-1.1219E-05
9.6591E-03
7.8096E-05
CUMULATIVE CHANSES
VOLUME IN <-)
VOLUME OUT
9.1921E+00
-1.9089E-03
D-13
-------
STORAGE 9.1798L+00
ERROR 1.7465E-02
RELATIVE ERROR 0.001903
EXECUTION ENDED NO*MALLY - 0 - WARNINGS
09/30/83 I0:0t:58 GCAGD.SOIL.LIB.DATA(GRIOl)
50 50.0
D-14
-------
09/30/83 10:05:17 GCAGD.SOIL.LIB.DATA(PROPl)
1
2
3
4
5
6
7
8
9
10
11
12
13
1*
15
16
17
18
19
20
21
22
23
2*
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
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.230E-05
1.230E-05
1.230EXJ5
1.230E-05
1.230E-05
1.230E-05
1.230E-05
1.230L-05
1.230E-05
1.230E-05
1.230E-05
1 .230C-C5
1.230E-05
1.230E-05
1.230E-05
1.230E-05
1.230E-05
1 .230E-05
1.230E-05
1.230E-05
1.230E-05
1.230E-05
1 .230E-05
1.230E-05
1.230E-05
1.230E-05
1 .230E-05
1.230E-05
1 .230c-05
1 .230E-05
1.230E-05
1.230E-05
1 .230E-05
1.230L-05
1.230E-05
1.230E-05
1.230E-05
1 .230E-05
1 .230E-05
1 .230t-05
1 .230E-05
1.230E-05
1.230E-05
1 .230L-D5
1.230E-05
1.230E-05
1 .230E-05
1.230E-05
1 .230E-05
1 .230E-05
0.4950
0.4950
0.4 950
0 . 4 95 C
0 . 4 95 C
0.495P
0.4950
0.4950
C . 4 95 0
0.4950
0 . 4 95 0
0.4950
0.49SD
0.4950
0 . 4 95 C
0 . 4 95 0
0 . 4 95 C
0.495P
0.4950
0.495C
0 . 4 95 0
0 . 4 95 0
u.4950
0.4950
0.4950
0.4550
0 . 4 95 0
0.4950.
0.4950
0.4950
0.4950
0 . 4 95 0
0 . 4 95 0
0 . 4 95 0
0.4950
0.4950
C.495C
0.4950
0 . 4 95 0
0.4950
0.4950
0.4950
0.4950
0 . 4 95 0
0.4950
0.4950
C.4950
0.4950
C.4950
C.4950
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1. 00
-1. 00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-l.CO
-1.00
-1.00
-1.00
-l.OC
-1.00
-l.CO
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
-l.DO
-1.00
-1.00
-1.00
-1.00
-1.00
-1. 00
-1.00
-1.00
-1.00
-1.00
-1.00
-1.00
D-15
-------
09/30/83 1C:05:21 GCA&D.SOIL.LIB.DATMIMTIAL3)
-600.0
09/30/83 10:05:29 GCA&D.S01L.LIP.OATA
-------
APPENDIX E
COMPUTER PROGRAM FOR GARDNER'S ANALYTICAL SOLUTION
E-l
-------
09/26/83
08:40:42
GCAGO.GARDNER.FORT.DATA
UUUU U J1U
OOODP020
00000^30
00000040
OOGOCG50
OOOOOCfaO
00000070
00000080
00000090
00000100
00000110
?OOOC120
00000130
00000140
9000015C
00000160
00000170
00000160
n n ft n n i on
UUUUU17U
00000200
00000210
00000220
00000230
00000240
00000250
0000026G
0000027C
00000280
00000290
00000300
00000310
00000220
00000330
00000340
00000350
00000360
00000370
00000380
00000390
00000400
00000410
00000420
00000430
00000440
00000450
00000460
00000470
C0000480
00000490
00000500
00000510
00000520
00000530
00000540
00000550
00000560
00000570
00000580
00000590
00000600
c
r
c
c
c
c
c
c
c
\_
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
GARDNER. FORT STtAQY STATE EVAPORATION IN UNSATURATEu SOIL
SEE J.R. GARDNER, SOM£ STEADY-STATE SOLUTIONS OF THE UNSATURATED
MOISTURE FLUb EMJATICN WITH APPLICATION TO EVAPORATION FROf,
A UATdR TABLE, SOIL SCIENCE, 85, 228-232, 1958.
STEADY STATE SOLUTION FOR SOIL UITH HYDRAULIC CONDUCTIVITY
FUNCTION AS FOLLOWS:
K = A S = - (PSD SUCTION
4
S + BETA
DAN GOODE JUNE 1983
IRO = 5
IPRT=6
IRSLT=9
REA3 SOIL FUNCTION PARAMETERS
READ(IRD,«) A,B,JT,FLUX
URITEIIPRT,2001> A,B,*T,FLUX
COMPUTE SOLUTION PARAMETERS
SGR2=1. 41421 2k
ALPHA=FLUX/A
BETA=ALPHA«b » 1.0
RO=(B£T A/ALPHA)** 3. 25
R02=RO*RO
R03=H02*RO
AINV=1.0/ALPHA
TERM1=1.0/(4.*R03*SQR2>
TERM?=1./(2.*R03»SQR2)
URITE(IPRT,2004) ALPH A, BETA ,R 0,R02 »R03,A INV ,TERM1,T ERM2
READ STEPPING FOR SOLUTION
READ PF,DPF»PFEND,MAX
URITt(IPRT,2002) PF,DPF,PFEND,MAX
PF=PF-DPF
DO 10 1 = 1, MAX
PF=PF+DPF
IF(PF.GT.PFEND) GOTO 999
S=10.0**PF
PSI=-S
S2=S*S
TERM3=RO*S*SQR2
ARGlr
-------
00000610
00000620
00000630
00000610
00000650
00000660
Q0000t70
000006&0
00000690
00000700
30000710
00000720
OC000730
0000074C
COOOG75C
COOOC160
OOOOC 770
00000760
00000790
00000800
00000010
00000620
00000o30
0000 0 640
00000850
0000086C
00000 870
00000 P80
1C
995
2C01
WRITE(IRSLT.2003)
CONTINUE
I«Z.PSItPF
STOP
FORf»ATF L'ND= • . OPF 1 0 . 2/
*• MAXIrtUM NUMBER UF p HINTS' f 5 (1H .). *MAX= • t 110///
*• I Z PSI PF S
•3 ARG1 AR32*«/>
2003 FCRMAT
------- |