-------
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» EXECUTION ENDED NORMALLY - «,I5,« - WARNINGS')
END
SUBROUTINE SPROP1 (PS I ,PSICR V, RMCUR V,RKCUR V,RMOIST,C ,RK ,
* RMSTL,CL,RKL,STAKK,SATK,POH,ISOIL,PSICRT,KODC>
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 2956.DO*(POR(L)-0.124DO)*PSILN**3/PSINEG/(739.DO»PSILN**4)**2
IF(JSOIL.EQ.2) CL(IP) -
* 6.379606*IPOR(L)-0.075DO)*PSINEG**2.96DO/
<(1.611D6+PSINEG*«3.9faDO)**2)
IF(JSOIL.EG.3> CL
C-6
-------
00003060
0000307C
00003080
oooosoyo
00003100
00003110
00003120
00003130
00003140
00003150
00003160
00003170
C0003180
000031*0
000032CO
00003210
00003220
00003230
00003240
C0003250
00003260
00003270
00003280
00003290
00003300
00003310
00003320
00003330
00003310
00003350
00003360
00003370
00003380
00003350
00003400
00003410
00003420
00003430
00003440
00003450
00003460
00003470
00003460
000034VO
00003500
00003510
30003520
20003530
00003540
00003550
00003560
00003570
00003580
00003590
00003600
nnnn^Ai n
uUUUJOXU
00003620
00003630
00003640
00003650
00003660
00003670
CBUG
6667
C
C
C
5
10
C
C
C
C
CDP
2C
C
C
C
CDP
30
C
C
C
CBUG
CBUG
CBUG
CBUG
CBUG
CBUG
CBUG
C
C
C
C
C* • • •
WRITE PS1NE5,PSILN,CLUP>.RKL.RMSTL(IP>
FORMA T< CP2F10.2,lF2£10.2fCPF10.4>
IF(K.EQ.l) GOTO 8
GOTO 1C
SATURATED NODE!
CONTINUE
RKSTL(IP)=POR(L>
CL< IP)=O.ODO
RKL(IP) = SATK (L)
If(K.EU.l) GOTO 8
CONTINUE
ENO OF LOOP
COMPUTE GEOMETRIC MEAN INTERBLOCK CONDUCTIVITY
12 = 0
DO 2D L=1.NUMEL
11=12+1
12=11+1
CHAMCic. NEXT CARD FOR DOUBLE/SINGLE PRECISION
STARK(L)=SURT(RKL(11 )*RKL(I2) >
CONTINUE
COMPUTE MEAN MOISTt K, AND C AT NODES
L2 = l
DO 30 I=2iMPMl
L1=L2+1
L2=L1+1
C(I)=(CL(L1)+CL(L2))/2.DO
RMOISTC I)=(RMSTL( L1)+RMSTL(L2))/2.DO
CHANGE NEXT CARD FOR DOUBLE/ SI NGLE PRECISION
RK(I)=SQRT(RKL(L1 )*RKL(L2)>
CONTINUE
USE ONE NODE VALUE AT TOP AND BOTTOM
C(1)=LL(1)
RMOIST(1)=RMSTL<1 )
RK«1) =RKLI 1)
C
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/« SPATIAL DISCRETIZATION P AR AMATERS • //
*' NUMBER OF NODE POINTS ' , 10 ( 1H. > , • 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»AT/»
*•
GARDNER OUTPUT*//
SOIL FUNCTION PARAMETERS'//
•• A»i 15Uh.) »»A = ».1PE10.2/
*' 6' t 15 (1H. ) t»B = » .1PE1C.2/
*• IE=TH TO WATER TAB LL'» t5 (1 H. ) , »»IT = • t OPF 1 0 .21
*• FLJX UPyARD».7( 1H. ) t»FLUX=» .1PL1C.2 )
2002 FOR«AT(//« SOLUTION STEPPING PAR AMLTERS* / /
*• BEGINNING SUCTICFJ' , 10(1H. ) . «PFr« , OPF10.2/
•• PF INCREMENT', 1 G(1H.) , • DPF= • , OPF1 0. 3/
*• FINAL PF«, 12<1H.) . • h>F 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
-------