EPA-600/2-77-1796
August 1977
Environmental Protection Technology Series
ION OF MINERAL QUALITY OF
IRRIGATION RETURN FLOW
Volume V. Detailed Return Flow
Salinity and Nutrient
Simulation Model
Robert S. Kerr Environmental Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Ada, Oklahoma 74820
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health-Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioeconomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the ENVIRONMENTAL PROTECTION TECH-
NOLOGY series. This series describes research performed to develop and dem-
onstrate instrumentation, equipment, and methodology to repair or prevent en-
vironmental degradation from point and non-point sources of pollution. This work
provides the new or improved technology required for the control and treatment
of pollution sources to meet environmental quality standards.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/2-77-179e
August 1977
PREDICTION OF MINERAL QUALITY OF
IRRIGATION RETURN FLOW
VOLUME V
DETAILED RETURN FLOW SALINITY AND
NUTRIENT SIMULATION MODEL
by
Marvin J. Shaffer
Richard W. Ribbens
Charles W. Huntley
Bureau of Reclamation
Denver, Colorado 80225
EPA-IAG-D4-0371
Project Officer
Arthur G. Hornsby
Source Management Branch
Robert S. Kerr Environmental Research Laboratory
Ada, Oklahoma 74820
ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OKLAHOMA 74820
-------
DISCLAIMER
This report has been reviewed by the Robert S. Kerr Environmental
Research Laboratory, U.S. Environmental Protection Agency, and
approved for publication. Approval does not signify that the con-
tents necessarily reflect the views and policies of the U.S.
Environmental Protection Agency, nor does mention of trade names
or commercial products constitute endorsement or recommendation
for use.
11
-------
FOREWORD
The Environmental Protection Agency was established to coordinate
administration of the major Federal programs designed to protect the
quality of our environment.
An important part of the Agency's effort involves the search for
information about environmental problems, management techniques and
new technologies through which optimum use of the Nation's land and
water resources can be assured and the threat pollution poses to the
welfare of the American people can be minimized.
EPA's Office of Research and Development conducts this search
through a nationwide network of research facilities.
As one of these facilities, the Robert S. Kerr Environmental
Research Laboratory is responsible for the management of programs to:
(a) investigate the nature, transport, fate and management of pollutants
in groundwater; (b) develop and demonstrate methods for treating waste-
waters with soil and other natural systems; (c) develop and demonstrate
pollution control technologies for irrigation return flows; (d) develop
and demonstrate pollution control technologies for animal production
wastes; (e) develop and demonstrate technologies to prevent, control
or abate pollution from the petroleum refining and petrochemical
industries; and (f) develop and demonstrate technologies to manage
pollution resulting from combinations of industrial wastewaters or
Industrial/municipal wastewaters.
This report contributes to the knowledge essential if the EPA is
to meet the requirements of environmental laws that it establish and
enforce pollution control standards which are reasonable, cost effective
and provide adequate protection for the American public.
William C. Galegar
Director
Robert S. Kerr Environmental
Research Laboratory
111
-------
PREFACE
This report is one of a set which documents the development
and verification of a digital computer modeling effort to predict
the mineral quality changes in return flows occurring as a result
of irrigating agricultural lands. The set consists of five separate
volumes under one general title as follows:
"Prediction of Mineral Quality of Irrigation Return Flow"
Volume I. Summary Report and Verification
Volume II. Vernal Field Study
Volume III. Simulation Model of Conjunctive Use and Water
Quality for a River Basin System
Volume IV. Data Analysis Utility Programs
Volume V. Detailed Return Flow Salinity and Nutrient
Simulation Model
This set of reports represents the culmination of an effort
started in May 1969 by an interagency agreement between the U.S.
Bureau of Reclamation and the Federal Water Pollution Control
Administration on a joint research proposal on the "Prediction
of Mineral Quality of Return Flow Water from Irrigated Land."
This research project has had three different project identifica-
tion numbers during the project period. These numbers (13030 EH,
EPA-IAG-048-(D), and EPA-IAG-D4-0371) are given to avoid confusion
on the part of individuals who have previously tried to acquire
project reports for the earlier project numbers.
IV
-------
ABSTRACT
A return flow quality simulation model is described which models the
plant-soil-aquifer system from the soil surface to a tile or open drain.
Processes simulated include evapotranspiration, unsaturated and satu-
rated water flow, solution-precipitation of slightly soluble salts, ion
exchange, ion pairing, nitrogen transformations, crop uptake of nitro-
gen, and the movement and redistribution of salts and nutrients.
The dynamic non-steady-state model predicts the concentrations of calcium,
magnesium, sodium, ammonium, bicarbonate, carbonate, chloride, sulfate,
N03-N, and Urea-N contained in soil, aquifer, and drain waters. Concen-
trations of organic-N; exchangeable calcium, magnesium, sodium, and
ammonium; and gypsum are predicted within the soil and aquifer materials.
Users' manuals for each basic subprogram are included, and a sample
problem illustrates the use of the model. Model output can serve as
input to the conjunctive use model described in Volume III; serve
as input to other models on a nodal or point source basis; or stand
alone depending on the type and scope of the particular study.
This report was submitted in fulfillment of project EPA-IAG-D4-0371
by the U.S. Bureau of Reclamation, Engineering and Research Center,
under the partial sponsorship of the Environmental Protection Agency.
Work was completed as of June 15, 1975.
-------
CONTENTS
Sections Page
Ifir TIT. -1 w... .1 IBI I!• Jl_ 11 • I , uBimiB-liLii ,,
I Conclusions 1
II Recommendations 2
III The Model 4
Introduction 4
Model Description - General 6
Model Description - Irrigation
Scheduling 6
Model^Description - Unsaturated
Flow 9
Model Description - Drainout 10
Model Description - Saturated Flow 11
Model Description - Interface
for Chemistry 12
Model Description - Unsaturated
Chemistry 12
Model Description - Saturated
Chemistry 13
Model Description - Drain Effluent
Prediction 14
IV Model Verification 16
V Model Strengths and Weaknesses 24
VI Summary of Required Input Data 39
Drainage Design 39
Saturated Material Properties 39
Unsaturated Material Properties 39
Soils Analysis 40
Crop Information 41
Water Applications 41
Fertilizer Information 41
Irrigation Water Analysis 42
Soils Temperature Data 42
Soils Coefficients 42
VII Interfacing with Conjunctive Use Model 43
VII
-------
CONTENTS - Continued
Sections
VIII User's Manual - Overall Model
Introduction
Subprogram Linkage
Bypass Options
Irrigation Scheduling
Unsaturated Flow and Unsaturated
Chemistry
Saturated Chemistry
Drain Effluent Prediction
Saturated Flow
Tips on Model Application
Use of Bypass Options
Intermediate Output
Printed Output Volume ,
Accuracy Versus Computer Run Times ..
Conversion to Other Computer Systems
IX User's Manual - Irrigation Scheduling
Program
Introduction
Data Input Card Sequence and Data
Limitations
X User's Manual - Unsaturated Flow Program .,
Introduction
Data Deck Structure
Node Spacing ,
Moisture Contents
Unsaturated Flow Properties
Initial Boundary Conditions ,
Consumptive Use Values ,
Depleted Moisture Option
Starting and Stopping Dates
Deficit Moisture
Mass Balance
44
44
44
44
44
45
46
46
46
47
47
54
54
55
55
56
56
56
70
70
72
80
80
81
83
84
85
86
86
87
Vlll
-------
CONTENTS - Continued
Sections
Page
Restart of Problem 90
Limitations 93
Output 93
Example of Subroutine PR0P 98
Campbell's Method for Determining
Unsaturated Flow Properties 103
XI User's Manual - Drainout Program 110
Introduction 110
Unit Numbers 110
Data Deck Structure Ill
Design Inputs 117
Selection of Drainage Case 119
Application of Method 122
Dynamic Equilibrium 123
Deep Percolation Data 123
Negative Starting Position , 124
Restart of Problem 126
Design Method 127
Input Checks and Error Messages 129
Output 129
XII User's Manual - Interface for
Chemistry Program 133
Introduction 133
Unit Numbers 134
Data Deck Structure 135
Tape Output Details 136
XIII User's Manual - Unsaturated
Chemistry Program 139
Introduction 139
Card Groups 139
Tape Inputs 153
Card Outputs 156
Tape Outputs 157
Hints on Program Use 158
Suggested Load Sequence 160
IX
-------
CONTENTS - Continued
Sections
XIV User's Manual - Saturated Flow Program ....
Introduction
Unit Numbers
Data Deck Structure
Standard Output
Input/Output Options
Grid Division and Geometry Modification
Limitations and Error Messages
Equations for Number of Nodes
Logical Error Messages
XV User's Manual - Saturated Chemisry
Program
Introduction
Unit Numbers
Data Deck Structure
Tape Input
Tape Output ,
Limitations
Output
Simulation of Unsaturated Zone
XVI User's Manual - Drain Effluent Prediction
Program
Introduction
Program Application
Quality of Leachate Water in Input
Data from Saturated Chemistry
Program is Input ,
Unit Numbers
Data Deck Structure
XVII References
XVIII Appendix
161
161
161
162
167
168
169
169
171
173
177
177
179
180
192
194
196
196
197
199
199
199
199
200
201
205
225
228
-------
LIST OF FIGURES
Number Page
1. Typical Vertical Section Through
Drainage F ield 7
2. Schematic Representation of Water
Quality Model 8
3. Diagram of Field Application 19
4. Predicted and Observed Water Quality -
Navajo Wash (Chemistry included) 20
5. Predicted and Observed Water Quality -
Navajo Wash (Chemistry omitted) 23
6. Model Bypass Configuration No. 1 48
7. Model Bypass Configuration No. 2 49
8. Model Bypass Configuration No. 3 50
9. Model Bypass Configuration No. 4 51
10. Model Bypass Configuration No. 5 52
11. Model Bypass Configuration No. 6 53
12. Modeled System - Unsaturated Flow 71
13. Data Deck Structure - Unsaturated Flow
Submodel 73
14. Punched Card Restart Deck - Unsaturated
Flow Submodel 91
15. Unsaturated Flow Properties - Exponential
Forms 100
16. Unsaturated Flow Properties 102
17. Data Deck Structure - Drainout Submodel 112
18. Program Criteria to Select Drainage Case 121
19. Organization - Daily Moisture Tape 125
20. Organization of Magnetic Tape Output -
Drainout Submodel 132
21. Card Groups in Unsaturated Chemistry
Submodel 140
22. Data Deck Structure - Saturated Flow
Submodel 163
23. Boundary Conditions - Saturated Flow
Submodel 165
24. Dimensionless Flow - Circular Tile Drains
Below Horizontal Water Table 170
25. Segmented Flow Tube in Saturated Region 178
26. Data Deck Structure - Saturated Chemistry
Submodel 181
27. Salt Inj ection from the Barrier 184
XI
-------
LIST OF FIGURES - Continued
Number
28. Organization of Tape 2 Output from
Unsaturated Chemistry Submodel 193
29. Organization of Tape 16 Output from
Saturated Chemistry Submodel 195
30. Data Deck Structure - Drain Effluent
Prediction Submodel 203
XII
-------
LIST OF TABLES
Number Page
1. Initial Conditions South Montezuma
Valley (Towaoc Area) 17
2. Average Annual Water Application Sequence 18
3. Present Observed and Predicted Conditions
(South Montezuma Area) 21
4. Present Observed and Predicted Conditions 22
5. General Model 25
6. Irrigation Scheduling Submodel 29
7. Chemical Equilibria Submodel 30
8. Nitrogen Transformation Submodel 31
9. Nitrogen Uptake Submodel 32
10. Salt Movement Submodel 33
11. Unsaturated Flow Submodel 34
12. Drainout Submodel 35
13. Saturated Flow Submodel 36
14. Saturated Chemistry Submodel 37
15. Drain Effluent Submodel 38
16. Drain Effluent Prediction Program 216
17. Drain Effluent Prediction Program 217
Xlll
-------
ACKNOWLEDGEMENTS
The foresight of Mr. John T. Maletic in supporting and promoting
modeling efforts from the earliest stages made a significant impact
on the rapid advances seen in this and other simulation models.
The pioneering efforts of Dr. Gordon R. Dutt in the development of
the chemical equilibria subprograms and his assistance with other
portions of the model are gratefully acknowledged.
Dr. Arthur W. Warrick and Mr. William J. Moore contributed the basic
portions of the Unsaturated Flow Program.
Messrs. Errol Jensen and Robert Julian of the Durango Planning Field
Division assisted in the model field verification study in the South
Montezuma Valley, Colorado area.
Messrs. Norman A. Roth, Roger W. Falkenstein, and Gary C. Shulstad
of the Missouri-Souris Projects Office developed the basic input
data used in the sample model runs contained in the Appendix.
The major sponsorship of this work was provided by the Mid-Pacific
Region and Water Quality Office of the U.S. Bureau of Reclamation.
xiv
-------
SECTION I
CONCLUSIONS
The simulation model described in this volume is a highly flexible,
sophisticated tool which can be used to predict the effects of initial
soil and aquifer conditions, cultural practices, input water quality,
and other factors on water and soil chemistries throughout the plant-
soil-aquifer system. Output can be used as input to other models
which consider the combined effects of several point sources on stream
quality in a river basin or subbasin.
The model also has many stand-alone applications such as predictions
involving the reclamation of saline and sodic soils, drain spacing
design, salt and nutrient pollution studies on soils and aquifers,
consumptive use predictions, and irrigation scheduling.
Potential model applications (those which would require additional
programing or modifications) include prediction of crop response to
salinity, optimization of nitrogen fertilization schedules, irrigation
scheduling (including unsaturated flow and salinity), and prediction
of additional quality parameters such as pesticides, phosphates, and
trace elements.
The flexibility built into the model allows the user to select the
degree of modeling sophistication needed for a particular application.
This can range from rough approximations to complete usage of the
model's capabilities.
-------
SECTION II
RECOMMENDATIONS
1. The total simulation capabilities of this model should be utilized
only when a maximal effort is indicated by the data and type of applica-
tion. For reconnaissance and survey applications, numerous submodel
bypass options are available to the user.
2. Although the primary outputs from the model are flow and quality
predictions associated with a tile or open drain, other model applica-
tions involving predictions within the plant-soil-aquifer system are
extremely useful.
3. This model (Volume V) should be interfaced with the conjunctive
use model (Volume III) when flow and quality predictions are wanted
for a river basin or subbasin (e.g. more than one node), when more
flow and chemistry detail is wanted in the unsaturated and saturated
zones, or when predictions of nitrogen transformations are desired.
4. Additional research is needed to develop and refine methodology
to select the input data used in this and other simulation models.
5. Additional field verification must be undertaken on this and
other models before they can be applied with a high degree of confi-
dence. The current confidence level of this model can be rated as
moderate with the accuracy of field predictions (areas on the order
of thousands of acres) appearing to be within about 10 to 15 percent
of observed data. Future verification studies could raise or lower
this figure.
6. The user should read and thoroughly understand the basic meth-
odology associated with the model before applying it.
7. An important application of this model should be its use as a
basis to develop more sophisticated simulation tools.
8. Although this model can (and has) been converted to other com-
puter systems, consideration should be given to utilization of the
Bureau of Reclamation CYBER 74-28 system via remote terminals.
9. This model represents a "first generation" attempt at modeling
the plant-soil-aquifer system. A "second generation" model can and
should be developed at an early date. The following features should
be incorporated into this more advanced modeling tool:
a. Unsaturated and saturated flow in three dimensions.
b. Multiple soil horizons in the flow submodels.
2
-------
c. A moving water table in the unsaturated flow submodel.
d. An improved inorganic chemistry submodel which incorporates
additional ion pairs and solid phases under certain conditions.
e. An improved denitrification submodel.
f. An improved plant uptake of nitrogen submodel including a
variable root distribution and improved uptake functions.
g. Additional submodeling to allow simulation of acid soil
chemical reactions.
h. Temperature effects in the inorganic chemistry submodel.
i. Generalization of certain reaction rate equations in
nitrogen transformation submodel.
j. Addition of submodels which simulate phosphorous, pesticide,
and trace element reactions and movement in soils.
k. Direct modeling of mechanical dispersion.
1. Fertilizer scheduling capability.
m. A submodel simulating crop response to salinity.
-------
SECTION III
THE MODEL
INTRODUCTION
During the past several years there has been a growing awareness
of the environmental problems created by man's present and past
activities. Prominent are those concerning water quality. When-
ever water is withdrawn from a natural supply and used or consumed,
its characteristics are changed. In particular, irrigation water
usually experiences significant alterations in passing through the
soil system before it is returned to the surface, either through
natural or artificial drainage pathways.
The nature of these alterations depends on many interrelated char-
acteristics of the soil-water system. Included are physical prop-
erties governing percolation of water through a porous media and
chemical characteristics which determine the nature of the chemical
pathways. Both are essentially beyond man's control. In contrast,
the design of artificial drainage systems, selection of crop types,
and fertilization and irrigation practices are largely within man's
control.
Viewed as a whole, the soil-water system is quite complex. In order
to understand the processes involved, workers have traditionally
dissected the system for closer inspection. Thus, the agronomist
is mainly concerned with crop productivity and factors affecting it
and is only indirectly concerned with movement of water beyond the
root zone. In contrast, the drainage engineer is primarily concerned
with estimates of deep percolation and the need to keep the water
table below the root zone. In reality, activities in both saturated
and unsaturated zones are related. Changes in cultural and irri-
gation practices affect the entire system to some degree, although
transformations in time and space may mask the cause-effect rela-
tionship. However, in order to evaluate these relationships, it is
essential that the individual components again be interfaced and
integrated.
When new or former dryland farms are brought under irrigation, it
is desirable to forecast the quality of the return flows. If sim-
ilar areas have already been developed, their performance can be
used to make this estimate. Where this is not possible, extrapola-
tion from dissimilar areas can lead to gross errors. As an alter-
native, a simulation model of the soil-water system, based on the
essential physical processes in it, can aid in predicting return
flow quality.
-------
The estimated results from a simulation of the system should not be
considered or represented as exact. There are three major sources
of uncertainty:
a. Model approximations and assumptions. These are made to
simplify the model. They may be dictated by a lack of basic
knowledge, by a lack of required data, or by the accuracy of
the results required.
b. Field variability. The parameters of a natural system may
possess significant spatial and temporal variations. Costs of
data collections generally prevent acquisition of sufficient infor-
mation to define it completely. Subjective estimates and intuitive
judgment are often employed to keep costs down.
c. Stochastic processes. For example, rainfall, solar radia-
tion, and temperature processes, are not completely understood
and apparently include random components. These affect plant
evapotranspiration and growth processes, as well as other por-
tions of the system.
Despite these uncertainties, a simulation model can be used to make
the most of the available data. It becomes another tool available to
the decisionmaker in assessing the impact of an irrigation project.
The detailed return flow quality simulation model described in this
volume was developed over a period of about 7 years and has had
several contributors (see acknowledgements). The model simulates
chemical and physical processes associated with agricultural lands
drained by subsurface tile drainage systems. The simulation begins
with field applications of water, salts, and nutrients (nitrogen)
and ends with predictions of flow and water quality from the drains.
In addition, the model can be applied to areas with surface drainage
systems and can provide predictions of quality and flow parameters
at various additional points within the plant-soil-aquifer system.
The background, theory and development of the various subprograms
within the model have been published elsewhere and need not be
repeated here.
Numerous references are included in the brief descriptions of each
subprogram. The reader is advised to thoroughly familiarize himself
with these materials before attempting to utilize the model.
Detailed user's manuals and a sample problem have been included to
facilitate application of the model.
-------
MODEL DESCRIPTION - GENERAL
Consider a typical vertical section through a drainage field as shown
in Figure 1. Water is applied at the land surface, either as irriga-
tion or precipitation. A portion evaporates or becomes surface run-
off. The remainder infiltrates the soil system and becomes available
for use by the plants or crops. In addition to evapotranspiration,
there is nitrogen uptake by the plants. As water percolates verti-
cally to the water table, chemical reactions occur. Upon reaching
the water table, flow takes place under saturated conditions, even-
tually reaching the drains. Chemistry changes continue to occur in
this saturated system.
This brief description leads to the schematic representation of the
water quality model as shown in Figure 2. Each program part or block
is identified with a major component of the physical system. The
conceptual model was translated to mathematical notation and programed
in Fortran IV for solution by a large, third-generation computer. All
programs are operational on the Control Data Corporation CYBER 74-28
computer located at the USER Engineering and Research Center in Denver,
Colorado. A brief description of each block follows.
MODEL DESCRIPTION - IRRIGATION SCHEDULING
The timing and amounts of irrigation are required for input to the
Unsaturated Flow Program. Deep percolation volumes are needed by
the Drainout Program. This information may be developed subjectively
based on past experience, by consideration of cultural practices in
similar areas, or by using manual or automatic irrigation scheduling
techniques. The latter is desirable where rainfall contributes sig-
nificantly in meeting plant requirements.
The philosophy involved in irrigation scheduling is to keep a daily
account of soil moisture in the root zone of the crop. This is accom-
plished by considering major inputs (irrigation and precipitation) and
withdrawals (crop evapotranspiration and deep percolation). At pres-
ent, a separate program is used. Since the simple bookkeeping pro-
cedures used in this program require estimates of deep percolation
losses and ignore possible contributions of water to the root zone
by capillary flow, it is anticipated that the scheduling techniques
will be incorporated in the Unsaturated Flow Program. This will
avoid duplication of elements common to both programs, avoid incom-
patibilities, and speed the scheduling process.
Crop evapotranspiration is computed using the Jensen-Haise method
(7, 8). Evapotranspiration is increased for 3 days following a rain-
fall or irrigation based on a procedure used by Jensen (8, p. 956).
This curve relates evapotranspiration at various growth stages to
-------
Land surface
Crops (evapotranspiration or consumptive use)
Bounding streamline
e>
LU
CO
Impermeable barrier
FIGURE I. TYPICAL VERTICAL SECTION THROUGH DRAINAGE FIELD
-------
DRAINOUT
IRRIGATION
SCHEDULING
UNSATURATED
FLOW
SATURATED
FLOW
INTERFACE
FOR
CHEMISTRY
UNSATURATED
CHEMISTRY
SATURATED
CHEMISTRY
DRAIN
EFFLUENT
PREDICTION
FIGURE 2. SCHEMATIC REPRESENTATION
OF WATER QUALITY MODEL
-------
the estimated potential evapotranspiration. If soil moisture is
depleted below the allowable depletion percentage, evapotranspiration
is reduced linearly, going to zero when all soil moisture is used.
In addition to irrigation scheduling information, average evapotran-
spiration rates on a semimonthly basis are computed for input to the
Unsaturated Flow Program. Rainfall also is accumulated on a semi-
monthly basis and applied on the 15th and last day of the month.
MODEL DESCRIPTION - UNSATURATED FLOW
The unsaturated flow submodel describes the infiltration, redistri-
bution drainage, and plant root extraction of soil moisture by a
growing crop (5, pp. 8-21). Darcy's Law is assumed to apply in
relating velocities to total hydraulic heads. Two modifications are
made:
a. The permeability (hydraulic conductivity) is taken as a func-
tion of moisture content. This accounts for the reduced area
available for flow at lower moisture contents.
b. The pressure (tension) head is taken as a function of mois-
ture content. This accounts for the fact that different forces
predominate at different moisture levels. Thus at high levels
gravity is most important while at low levels capillary or molec-
ular forces have greatest effect.
Hysteresis in both the tension head and conductivity functions is
ignored. Multiphase flow (simultaneous movement of air, water vapor,
and water) is not treated. Flow is assumed to be independent of
quality, although the chemical processes depend on the flow.
Since both the conductivity and head are highly nonlinear with mois-
ture content, solution of the flow equation is accomplished using
numerical methods. An implicit finite difference scheme is employed
and flow is treated as one dimensional in the vertical direction.
Water applications are applied on the proper date and represent "effec-
tive" values. Thus, surface runoff and evaporation from depression
storage are subtracted from the total amounts. Withdrawal of water by
the plants is accomplished on a macroscopic basis using a "sink" term
in the governing flow equation. Thus, water is extracted continuously
from the land surface to the bottom of the root zone in proportion to
the root distribution. Details of flow to individual roots are not
considered. Plant evapotranspiration is estimated by the Jensen-
Haise method utilized in the Irrigation Scheduling Program or by using
the Blaney-Criddle method. At present, root distributions are taken
as constant with time and are chosen to reflect the mature growth stage.
-------
Water to meet the specified evapotranspiration is removed node by node
as long as moisture levels are above a specific lower limit, usually
taken as the wilting point. If levels fall below this limit at a node
no water will be withdrawn from that node and the evapotranspiration is
reduced. The reduction is accumulated and listed as the "deficit mois-
ture." No consideration is made for effects due to permanent wilting.
When water contents again go above the wilting point, the plants use
water at the specified evapotranspiration rates. Effects of salinity
on evapotranspiration are also ignored.
The lower boundary of the one-dimensional column is taken as the water
table, which for purposes of the Unsaturated Flow Program is taken as
the mean annual location of the water table. Hydraulically it is
treated as a zero pressure surface maintained at the saturation mois-
ture content. The upper boundary is treated as a specified head bound-
ary when water stands on the surface after an irrigation and as a zero
flow boundary when all water from a surface application has infiltrated.
Moisture contents and water movements are then used as inputs, via
the interface program, to the Unsaturated Chemistry Program. Water
crossing the lower or water table boundary also acts as input to the
Drainout and Saturated Flow Programs.
MODEL DESCRIPTION - DRAINOUT
The quantity of water crossing the water table during a given day, as
computed by the Unsaturated Flow Program, is used to calculate the
position of the water table and the drain discharge as a function of
time. Standard Bureau methods for drainage system design are used.
These are described by Dumm (1,2,3).
Basically, the water table shape is chosen as a fourth-degree parab-
ola. Daily increments of flow across the water table are used in
conjunction with the specific yield to compute instantaneous rises
(if water moves from the unsaturated to saturated system) or instan-
taneous declines (if water moves from the saturated to unsaturated
system by capillary flow). Equations describing the flow of water
through the saturated system to parallel circular drains are then
used to compute drain discharge and the fall of the water table due
to drainout during the day.
Since the formulations are based on the Dupuit-Forchheimer assumptions
which treat flow as parallel and neglect the effects of convergence,
corrections for the depth of flow are made. Hooghoudt's equivalent
depth is employed for this purpose. Approximate equations developed
by Moody(ll) are used to simplify the program.
Discharge values, computed on a daily basis, are accumulated to yield
monthly and yearly values. By considering initial and final water
10
-------
table positions for each year, changes in storage during the year may
be computed. Since the amount of water crossing the water table is
known, application of a mass balance equation permits a semi-independent
determination of the annual drain discharge. If this differs signifi-
cantly from that computed using the drainout method, all daily discharge
values are adjusted to yield mass balance. This approach is used
because the discharge formulas are based on the first few terms of an
infinite series and the constants are increased for design purposes
(6, Eq 54, p. 41). Since the adjustment coefficients may differ from
year to year, daily discharge values may exhibit discontinuities
between years.
The resultant adjusted monthly discharge values are then used as input
to the combining or drain effluent programs.
MODEL DESCRIPTION - SATURATED FLOW
The total amount of deep percolation computed by the unsaturated flow
program is converted to an average annual rate for input to the satu-
rated flow program. Using the geometry of the drainage system, a
steady-state potential theory solution is employed to define stream
paths or lines to the drain (12, Section 2). These solutions remove
the limitations imposed by the Dupuit-Forchheimer assumptions.
A number of stream tubes may be defined to represent the two-
dimensional flow system. Since an average annual seepage rate is
used, a mean or average stream tube system results.
Because the equations defining the potential and stream functions are
in terms of infinite series involving trigonometric, hyperbolic, and
exponential functions, a node system is superimposed over the drain-
age system. Stream functions are computed at the nodes and required
values found by linear interpolation. The corresponding potentials
are then computed. With streamlines located at the hydraulic center
of each tube and the potentials known, Darcy's Law is employed to
calculate velocities of flow between adjacent points on the stream-
line. Using the known distances and porosities, these velocities
are converted to travel times. Integrating or summing the travel
times then yields the total time required for water crossing the water
table to arrive at the drains assuming piston displacement.
The water table position is also computed and the total volume of aqui-
fer contributing to one drain is found by numerical integration. This
volume is then distributed between individual tubes using the computed
travel times as the basis.
The pore volumes are used for input to the saturated chemistry pro-
grams and are necessary in calculating the soil mass and associated
11
-------
pore volumes in the control volumes used for calculating chemical
changes. Each stream tube is treated independently in the saturated
chemistry program and the results combined in the drain effluent pro-
gram. This combining is based on the relative travel times for piston
displacement in each tube.
MODEL DESCRIPTION - INTERFACE FOR CHEMISTRY
Moisture contents and movements from the unsaturated flow program
are computed at the nodes used in the finite difference equations.
In contrast, the unsaturated chemistry programs use a soil segment
concept in which each segment is assumed to possess homogeneous char-
acteristics. Because the segment centers and nodes of the unsat-
urated system may not coincide or be identical in number, an inter-
face program is necessary. This program computes average moisture
contents in each soil segment as well as movements of water between
soil segments for a given time interval.
MODEL DESCRIPTION - UNSATURATED CHEMISTRY
The Unsaturated Chemistry Program simulates a number of important
biological, chemical, and physical processes (5, pp. 21-79). These
include:
a. Nitrogen transformations, including the hydrolysis of urea,
mineralization-immobilization of organic-N and ammonium, nitri-
fication of ammonium, and immobilization of nitrate. Because
the rates of these transformations are slow (on the order of days
or weeks), a kinetic approach is used.
b. Salt or inorganic chemistry, including ion exchange, solution-
precipitation of slightly soluble salts, the formation of undis-
sociated ion pairs, and the bicarbonate buffer system. Since
reaction times are rapid (on the order of seconds or minutes),
chemical equilibrium conditions are assumed.
c. Unsaturated movement and redistribution of soluble constituents.
d. Crop uptake of nitrogen.
The nitrogen rate equations were developed by applying system analysis
concepts to define and limit the system and to establish pertinent
variables and reaction pathways. Equations based on these variables
were derived using multiple-regression analysis, reaction rate theory,
and statistical thermodynamics(5,13,14). Variables include ammonium,
nitrate, organic-N, and urea concentrations; temperatures; carbon-
nitrogen ratios; ionic strength; and soil moisture contents.
12
-------
Chemical reactions included in the salt (inorganic) chemistry portion
are the dissolution and precipitation of gypsum(5,16,17,21) and calcium
carbonate(5,21) undissociated ion pair reactions for calcium sulfate
and magnesium sulfate(5,16,17,21), calcium-magnesium ion exchange(5,17),
calcium-sodium exchange(5,21), sodium-ammonium exchange(S), and
pCOa - Ca++ - HCOa interactions(5). Constituents presently considered
in this portion of the model are nitrate, ammonium, calcium, sodium,
magnesium, bicarbonate, chloride, carbonate, and sulfate.
Solutes contained in the soil water are assumed to move with the
water into and/or from adjacent segments (5, p. 23). The volume of
water which moves and the mean water content of each segment (cell)
over a time step are inputed from the Unsaturated Flow Model via the
Interface Program. Because each cell is spatially homogeneous,
solutes transferred into a segment or cell are mixed with those
already in the cell.* This transfer may be either in an upward or
downward direction. After each mixing operation, the soil solution
phase is equilibrated with the solid and exchangeable phases. Since
a number of cells, usually 10 or more, are used, this numerical tech-
nique tends to simulate physical or mechanical dispersion. This type
of dispersion arises from tortuosity of the flow paths whereby flow
splits and branches in going around individual soil particles. Molec-
ular diffusion resulting from molecular energy movements is ignored
and is generally significant only when little or no water movement
occurs.
MODEL DESCRIPTION - SATURATED CHEMISTRY
In addition to the quantity and quality of leachate water predicted
by the Unsaturated Chemistry Program, the Saturated Chemistry Program
accepts the stream tube volumes determined in the Saturated Flow Pro-
gram as input. Each tube is treated independently and is subdivided
into a number of segments of equal volume. Leachate water is accu-
mulated until it equals the pore volume in a single segment. It is
then assumed to move by piston displacement through successive seg-
ments until it reaches the drain. (4) After each displacement, the
solution phase is equilibrated with the solid and exchangeable phases.
Lateral mixing between tubes and longitudinal mixing in adjacent seg-
ments is ignored. Thus, lateral dispersion and diffusion are ignored.
For saturated systems this is considered realistic.
* This technique applies only to the unsaturated zone. A different
procedure is used in the saturated region.
13
-------
Chemical transformations are computed after the slug of water enters
each segment. They are identical to those in the salt chemistry por-
tion of the Unsaturated Chemistry Program except that denitrification
is simulated as the deep percolation water crosses into the saturated
zone. Denitrification is assumed to occur at a rate which is temper-
ature dependent but independent of substrate (i.e., N03) concentrations
or zero order above a "saturation" level. This level has not been
well defined but probably is on the order of a fraction to a few ppm
N03-N. Below this level, the reaction rate becomes first order and
decreases as the N03-N concentration diminishes.
Initially, the entire saturated zone is considered as homogeneous.*
Laboratory analyses of soil samples in this zone are used to define
initial conditions. If more than one sample is available, analyses
may be averaged or a single typical analysis used.
As each slug of water enters the segment adjacent to the water table,
the slug contained in the segment adjacent to the drain is considered
totally displaced. The quality of water discharging from the given
stream tube is now that of the water in the segment adjacent to the
drain after all reactions have occurred. It remains at this quality
until another slug enters the tube at the water table. Individual
tubes are combined in the Drain Effluent Prediction Program.
Although the saturated system is initially homogeneous, the passage
of water down individual tubes soon produces heterogeneous conditions.
As each slug moves from segment to segment, its quality changes. The
chemistry of the soil complex also undergoes corresponding changes.
As more tubes and segments are used, the continuous two-dimensional
system is more closely approximated.
MODEL DESCRIPTION - DRAIN EFFLUENT PREDICTION
The Drain Effluent Prediction Program accepts as input the quality
of each slug discharging into the drain from individual tubes, the
mean travel times for piston displacement, and the monthly volume
of drain discharge. Travel times are used to calculate arrival times
of each slug at the drain on a tube-by-tube basis. Using these times,
the resultant quality of the drain effluent is obtained by combining
the appropriate slugs from individual tubes. This gives a varying
* The program will accept a different soil analysis for each segment.
However, because stream tubes curve and segments are of different
lengths, a considerable effort is required to use analyses by depth.
It is anticipated that future program additions will automate this
task.
14
-------
drain effluent quality under steady discharge. The monthly volumes
are then used to convert these results to the varying discharge pat-
tern with the corresponding monthly qualities. The accumulated salt
load out of the drain and the rate of removal are also computed.
All results are given for individual constituents and for the total
load. The total load is taken as the sum of the individual parameters
and approximates the total dissolved solids. Constituents presently
considered are nitrate-nitrogen, ammonia-nitrogen, calcium, sodium,
magnesium, bicarbonate, chloride, carbonate, and sulfate. Output
may include cathode ray tube plots, printed listings, and punched
cards. Results are on a "per-acre" basis. The assumption is made
that no chemical reactions occur after the mixing process in the
drain.
15
-------
SECTION IV
MODEL VERIFICATION
Verifications of the various subprograms within the model are included
in the literature pertaining to each program. Appropriate sections
describing each subprogram should be consulted for specific citations,
Verification of the model as applied to the field situation has been
obtained for an irrigated area in southwestern Colorado. An area of
about 2,200 acres has been irrigated for about 75 years. Soil samples
from an adjacent unirrigated area were utilized to estimate the ini-
tial soil conditions. The soil analyses from two sites were weighted
with depth to obtain mean soils data. These data appear in Table 1.
Profile B3 represents soils closest to the associated stream (Navajo
Wash) while profile B2 represents conditions farther from the stream.
Selection of sampling sites was based entirely on the experience of
field personnel.
Model bypass configuration No. 1 (see Section VIII, Figure 6) was
selected with the Irrigation Scheduling Program bypassed. An annual
average irrigation and precipitation input sequence was derived for
the simulation period (Table 2). This sequence was applied to the
irrigated area for the 75-year period.
A typical cross section through the irrigated area is shown in
Figure 3. The diagram includes the five stream tubes used in the
simulation together with mean travel times for each. Note the rela-
tively high travel times for flow through this aquifer.
Predicted and observed values for constituent and TDS concentrations
with time are given in Figure 4. Predicted and observed mean soil
profile data after about 75 years of irrigation appear in Table 3.
These soil profiles were selected in approximately the same relative
locations as the input profiles. The agreement between prediction
and observation was on the average within about 8 percent for soil
water quality and about 12 percent for stream quality. More detail
in the simulation such as additional soil profiles, layering with
depth, unsaturated flow and chemistry, additional ion pairs, and
better estimates of C02 partial pressure would probably improve the
agreement.
There has been some difference of opinion and confusion concerning
the importance of including the chemical reactions in field simula-
tions. (23, 24) To test the significance of the chemistry in this
study, computer runs were repeated in all respects except that the
subroutines which simulate the chemistry were bypassed. In effect,
16
-------
Table 1. INITIAL CONDITIONS
South Montezuma Valley
(Towaoc Area)
Soil profile B2 B3
Ca++
Na+
Mg++
HC03~
C03~
Cl"
80,,=
SUMS
24
52
57
1
0
8
118
260
8,369
29
110
51
2
0
32
145
369
11,940
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Mg/L
17
-------
Table 2. AVERAGE ANNUAL WATER APPLICATION SEQUENCE
Application
No.
Depth
meq/L
(Cm)
Ca++
Mg**
Na+
HC03-
C03"
S0^s
ci-
N03~
1
1.37
0.75
0.09
0.04
0.39
0.00
0.60
0.01
0.00
2
0.33
13.26
9.68
1.86
10.49
0.00
7.18
1.44
0.05
3
1.60
22.77
5.54
4.13
20.42
0.00
9.48
3.47
0.08
4
1.14
23.85
5.91
4.41
21.41
0.00
9.76
3.75
0.09
5
1.02
23.85
5.13
4.23
20.75
0.00
9.58
3.57
0.08
6
1.27
20.84
4.98
3.66
17.75
0.00
9.01
3.09
0.08
7
0.51
19.34
4.51
3.29
16.53
0.00
8.63
2.82
0.07
8
0.38
19.24
4.60
3.38
16.90
0.00
8.73
2.82
0.07
9
0.13
13.05
2.63
1.78
9.95
0.00
7.14
1.41
0.05
-------
1800 Ft.-
-Land surface
Mean travel time
59 Yrs.
FIGURE 3. DIAGRAM OF FIELD APPLICATION
-------
1000 -
(ppm)
1000 -
900 -
800 -
700 -
600 -
500 -
400 -
300 -
200 -
100 -
IDS •
S0 = •
Mg-n- x
Predicted-x-
10 20
30 40 50
TIME (YEARS)
60
Na+
o
TIME (YEARS)
HCOj a
10
20
30
40
i
50
60
70
I
80
TIME (YEARS)
FIGURE 4. PREDICTED AND OBSERVED WATER QUALITY
NAVAJO WASH (CHEMISTRY INCLUDED)
20
-------
Table 3. PRESENT OBSERVED AND PREDICTED CONDITIONS
(South Montezuma Area)
M2 (near Wash) M3 (farther from Wash)
Soil profile Observed Predicted Observed Predicted
Ca++ 25
Na+* 50
Mg"1"1" 29
HCOa" 2
co 3- o
Cl~ 3
S0k~ 97
206
SUMS
6,885
28
27
43
3
0
4
89
194
6,298
23
61
59
2
0
10
131
286
9,342
27
49
60
2
0
7
123
268
8,668
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Mg/L
Percent error 8.5 7.2
Mean
percent error 7.9
the applied waters were concentrated according to the consumptive
use and routed through the system. Initial conditions were identical
to those used previously.
Predicted and observed soil profile data after 75 years are given in
Table 4. The prediction error increased from 8.5 to 48.2 percent in
the soil closer to Navajo Wash and from 7.2 to 20.9 percent at the
more distant soil location. Since the soils closest to the drain are
more highly leached, the conclusion can be made that the chemical reac-
tions become increasingly more important as leaching progresses.
However, even in the soils more distant from the discharge point,
the improvement in predictive capability would easily demand inclusion
of the chemistry in the simulation.
21
-------
Table 4. PRESENT OBSERVED AND PREDICTED CONDITIONS
South Montezuma Area
(Chemical Reactions Bypassed)
M2 (near wash)
Soil profile Observed Predicted*
M5 (farther from wash)
Observed Predicted*
CA-1"1"
Na+
Mg++
HC03
CO 3=
Cl
so»r
SUMS
25
50
29
2
0
3
97
206
6,885
20
15
17
13
0
4
35
104
3,566
23
61
59
2
0
10
131
286
9,342
25
43
48
6
0
7
98
227
7,388
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Meq/L
Mg/L
20.9
Mean
percent error
34.6
Predicted and observed data for water quality in Navajo Wash with the
chemistry omitted from the calculations appear in Figure 5. These
results indicate the prediction error increased from about 12 percent
without the chemistry considered to about 28 percent without the chem-
istry after 70-75 years of irrigation.
A comparison of the predicted curves in Figure 5 with those in Figure 4
shows the Figure 5 TDS predictions are initially, high across the
Figure 4 TDS curve, and are low in the 70-75 year period. This empha-
sizes the variable conclusions which can be drawn concerning the impor-
tance of the chemical reactions. It appears from this analysis that
the chemistry is extremely important in most cases and should be
included in the simulation.
22
-------
(ppm)
4000 -i
3000 H
2000 H
1000
-IDS
10
20
TIME (YEARS)
30 40 50
TIME (YEARS)
TDS_«
so;«
Mg++X
No+A
cro
70
80
(ppm)
1000 -
900 -
800 -
700 -
600 -
500 -
400 -
300 -
200 -
100 -
IO
—I 1 1 1—
20 30 40 50
TIME (YEARS)
60
70
80
FIGURE 5. PREDICTED AND OBSERVED WATER QUALITY
NAVAJO WASH (CHEMISTRY OMITTED)
23
-------
SECTION V
MODEL STRENGTHS AND WEAKNESSES
The model contains numerous strong and weak points which may not be
readily apparent from the model descriptions and supporting docu-
mentation. The following set of tables contains 2 summary of major
strengths and weaknesses associated with the model. In addition,
comments are provided which should help clarify certain points. The
summary should not be considered an exhaustive treatise on the subject,
but merely an aid to better understanding the basic capabilities and
limitations of the model.
24
-------
Table 5. GENERAL MODEL
Strengths
Weaknesses
Comments
The model is not directly
tied to a specific irriga-
tion system and can be used
for almost any system by
using the programs properly
and interpreting the data
inputs properly.
The model must be applied
to each specific situation
encountered. Theoretically,
any change in inputs re-
quires a separate analysis.
Model is built on a modular
basis using subroutines and
the Fortran IV programing
language. It is easily
expanded.
Different systems such as
furrow and drip irrigation
must be approximated and
the assumptions "lived"
with.
Extrapolation to large areas
involves the problem of field
variability. It may be uneco-
nomical to collect the neces-
sary field data and make
model runs.
There are already a large
number of subroutines and
programs in the model.
Knowledge of programing is
also required.
The model is relatively
complicated and expensive to
run. Data requirements are
considered excessive by some
individuals.
Specifics can be added to
remove limitations and assump-
tions. Basic irrigation system
design is tied into the model by
efficiencies which are input to
the irrigation scheduling sub-
model. Efficiency also relates
cultural practices, crop and
soil types, climatic conditions,
and land slopes, etc.
The effect of areal variability
and the sensitivity of model
results to it is an item needing
further research. Practically,
results are more sensitive to
some inputs than to others. How-
ever, relative sensitivity can
change from problem to problem.
Special features of Fortran com-
pilers that are only available
at one computer center have gen-
erally been avoided.
Efforts to include items not
presently included will tend to
further complicate the model.
-------
Table 5. (continued) GENERAL MODEL
Strengths
Weaknesses
Comments
to
Ox
The model simulates many
important physical and
chemical processes which
have been identified
over the years and
expressed in equation
form.
The following quality para-
meters are included in the
model: Ca++, Mg++, Na"1",
NH.
NO,
co3, ci,
SO^, urea, organic-N,
exchangable Ca++ > Mg , Na+,
and NHt^"1", and the ion pairs
and MgSOt,.
The model has several appli-
cations other than predict-
ing the quality of return
flows. Those include soil
reclamation, irrigation
scheduling, drain spacing
design, and generating
inputs to other models.
The model is designed to
allow many computer runs to
be made in a short period of
time.
The model assumes away many
real life complications and
situations, making it unac-
ceptable to some theoreti-
cians. It is already so com-
plex as to be unacceptable to
many practitioners.
Additional quality parameters
such as phosphates, pesticides
trace elements, etc., are not
included.
Additional applications exist
which are still beyond the
immediate capabilities of the
model.
Few guidelines are provided
the user as to selection of
representative input data.
The model is in a "no-mans" land,
and may be unacceptable to the
two groups for opposing reasons.
A few individuals appear to be
suffering from future shock and
refuse to examine this and other
simulation models let alone
accept them.
The model is open ended and addi-
tional parameters can be added.
Other parameters such as ESP,
SAR, hardness, etc., can be cal-
culated from the existing
parameters.
The model will be improved and
expanded in the future.
Research is badly needed in this
area.
-------
Table 5. (continued) GENERAL MODEL
Strengths
Weaknesses
Comments
The model simulation begins
at the soil surface and ter-
minates at the drain (or
stream}.
Crop uptake of nitrogen is
simulated.
Chemical reactions are
taken as a function of
waterflow.
The model includes many
bypass configuration
options.
Model verification has been
obtained for several
applications.
The model utilizes chemical
equilibrium approach for
inorganic (salt) reactions.
Nitrogen transformations
are simulated utilizing a
kinetic approach.
Multiple drains/streams cannot
be mixed or combined in this
model.
A fixed crop root distribu-
tion is assumed. (One can be
entered for each crop.)
Uptake function is crude.
Waterflow is not considered
to be a function of the
chemistries.
The user may be confused as
to which option(s) to use.
Additional field verification
would be useful.
Some inorganic (salt) species
may not be in chemical
equilibrium.
Some rate equations are sta-
tistical and may not be valid
outside derivation data set.
Other models can accomplish this
function (e.g. - see Volume III)
A function(s) could be added
making the distribution variable.
An iterative approach could be
used to add this feature.
Additional model sensitivity
studies could help alleviate this
problem.
Adequate field verification will
require many years.
Model verification indicates this
assumption is excellent for the
species considered.
The transition-state nitrifica-
tion equation is general.
-------
Table 5. (continued) GENERAL MODEL
Strengths Weaknesses Comments
Hydrodynamic dispersion is Dispersion is not modeled Verification data indicate the
simulated using a numerical directly. model adequately simulates hydro-
scheme, dynamic dispersion.
Multiple stream tubes are The flow lines are taken as Uncertainties about aquifer prop-
considered in the saturated time invariant. erties make the current approach
zone. reasonable.
to
OO
-------
Table 6. IRRIGATION SCHEDULING SUBMODEL
Strengths
Weaknesses
Comments
Consumptive crop use is com-
puted using the Jensen-Haise
method.
The model makes use of con-
cepts long accepted by soil-
scientists such as field
capacity, wilting point,
and available moisture.
Considerable data are needed
to apply this method. Solar
radiation data is not col-
lected at all weather
stations.
The effect of soil moisture
on plant growth is ignored.
Indiscriminate use without
realizing the actual physical
processes involved can lead
to poor results in unusual
situations such as high water
table conditions.
The unsaturated flow model accepts
consumptive use computed by any
method. It could easily be
altered to accept basic data for
other methods such as Blaney-
Criddle. This is also true of
the scheduling model.
A plant growth model could be
added, provided sufficient data
are available to develop and use
it.
These methods have proven their
usefulness. The unsaturated flow
submodel tends to tie the physical
situation to the older concepts
of soil "constants."
-------
Table 7. CHEMICAL EQUILIBRIA SUBMODEL (Unsaturated Zone)
Strengths
Weaknesses
Comments
Chemical equilibria subpro-
gram simultaneously con-
siders several chemical proc-
esses important in many
Western and other soils,
including solution-
precipitation of lime and
gypsum, cation exchange, ion
pairing of CaSO^ and MgSO^,
and pC02-Ca++-HC03~
interactions.
pCOa is either calculated
as a function of soil mois-
ture content or entered into
the model as part of the
input data.
Activity coefficients are
computed using Debye-
Hiickel equation generalized
for monovalent and divalent
ions.
Chemical equilibria subpro-
gram ignores (1) other solid
phases such as Na2SOtt and
MgC03, (2) the formation of
additional ion pairs, and
(3) pC02-Ca+*-HC03~ interac-
tions in noncalcareous soils.
Ion exchange is limited to
Ca-Na, Ca-Mg, and Na-NH,/
exchange.
Other functional relation-
ships may be more useful
in computing pC02 under
situations such as within
a crop root zone.
Equation could be modified for
each specific ion. Uncharged
species are assumed to have
activity coefficients equal
to unity.
The model is not applicable where
solid phases other than lime and
gypsum are present at chemical
equilibrium and where additional
ion pairs become significant.
The model tests for supersatura-
tion with respect to lime, but
does not react C02 with the non-
calcareous system. The model does
not consider ion exchange in acid
soils.
The user must input pC02 values
within the root zone.
These approximations probably are
adequate in light of other
uncertainties.
-------
Table 8. NITROGEN TRANSFORMATION SUBMODEL (Unsaturated Zone)
Strengths
Weaknesses
Comments
Simultaneously considers
nitrification, urea hydrol-
ysis, mineralization, immo-
bilization, NH^* ion
exchange, and denitrifica-
tion (saturated zone only).
Nitrogen transformation rate
equations include important
variables such as substrate
concentrations, temperature,
soil moisture content, etc.
Nitrification is simulated
using either a statistical
equation, or a transition
state equation.
Transition state equation
includes pH, partial pres-
sure of 02> osmotic effects,
and NHtf* concentrations as
variables.
Ignores nitrogen fixation,
ammonia volatilization,
denitrification (unsaturated
zone).
Most rate equations were
derived statistically based
on experimental data from a
limited number of soils.
Extrapolation to other soils
may be dangerous.
Statistical equation may be
limited to certain soils.
Transition state equation is
general but applies only to
soils with fixed microbial
populations.
Microbial population is
assumed to be fixed in size.
Fixed nitrogen can be added to
soil as input data.
Transition state nitrification
equation is general.
Transition state equation could
be extended to include variable
microbial population size.
-------
Table 9. NITROGEN UPTAKE SUBMODEL
Strengths
Weaknesses
Comments
Nitrogen uptake by the crop
is allowed in the form of
NO," and NH, *.
Total nitrogen uptake can
either be inputed to the
model on a semimonthly
basis or computed by
assuming that N uptake is
proportional to consumptive
use. (The user supplies
the proportionality
constant.)
A fixed ratio between the
uptake of these two species
is used. Uptake of other
N species is not considered.
Additional nitrogen uptake
functions or schemes are not
considered.
A variable ratio could be added
at a future time. The major
N species taken as by plants are
N03" and NH,/.
The model is not yet particularly
sophisticated in its handling
of crop N uptake. This submodel
will be improved in the future.
-------
Table 10. SALT MOVEMENT SUBMODEL (Unsaturated Zone)
Strengths
Weaknesses
Comments
CO
GO
Unsaturated waterflow
volumes as computed by the
Unsaturated Flow Program
are used to move soluble
salts into and from each
soil segment (both in
upward and downward
direction).
Soil chemical reactions
are included in the leach-
ing process at each step
in time and space.
Mechanical dispersion is not
simulated directly.
Lateral movement of salts
is not considered.
Some waterflow rates may
be too high to allow chemical
equilibrium as assumed in the
model.
The use of several soil segments
per profile tends to introduce
numerical dispersion. This tech-
nique has yielded good correla-
tions with observed data in
unsaturated lysimeter tests.
-------
Table 11. UNSATURATED FLOW SUBMODEL
Strengths
Weaknesses
Comments
Unsatured waterflow is
simulated in the upward and
downward directions.
The bottom boundary is
taken as the water table.
Consumptive use is computed
using either the Jensen-
Haise or Blaney-Criddle
methods.
Hysteresis in the conduc-
tivity and diffusivity
relations is ignored, sim-
plifying program use and
reducing input requirements.
Two- and three-dimensional
flows are not considered.
Consequently, the model is
not directly applicable to
trickle and furrow irrigation.
The geometric location of the
water table is not changed.
Other methods may be more
desirable.
Hysteresis, especially in the
diffusivity, has been consid-
ered by many investigators as
important in "wetting-drying"
systems.
Additional dimensions could be
added but with increased computer
storage requirements and costs.
This tends to overestimate the
upward movement of water by
capillary forces when plant-
induced gradients produce such
flow. A moving water table loca-
tion could be modeled.
They can be added to the irriga-
tion scheduling or unsaturated
flow models or to both.
Inclusion in the model is pos-
sible, although it would compli-
cate programing, increase computer
storage, and lengthen runs. Ade-
quate field data to run the model
would be almost impossible to
obtain.
Water is extracted from the
root zone according to a
specified root pattern,
using a root extraction
term.
The root distribution is
constant with time.
A variable root pattern could
easily be added. It could be a
function of time or be based on a
more complicated model of plant
growth (as influenced by tempera-
ture, availability of water, and
the salinity of the soil water).
Field data would be costly.
-------
Table 12. DRAINOUT SUBMODEL
Strengths
Weaknesses
Comments
USER drainage equations are
utilized.
Submodel will either pre-
dict the response of a
given drainage system to
deep percolation inputs or
estimate the drain spacing
based on system physical
properties and desired
drainage depth.
Submodel computes the
monthly drain discharges.
Homogeneous, isotropic
conditions are modeled.
This allows use of rela-
tively simple mathematics
that are economical for
computer calculation.
Approximations associated with
USER techniques are included.
User may want discharges over
some alternate time period.
Real situation involves
heterogeneous conditions.
See comments.
USBR drainage equations have seen
a wide application in the United
States and foreign countries.
Model duplicates USBR method for
designing lateral spacings and
consequently is accepted by drain-
age personnel.
Daily discharges are actually
computed by the model and can be
listed or plotted as output. They
are available in the special plot
file generated by the program.
Weighted permeabilities with depth
are used to account for variabil-
ity of soils. The storage coeffi-
cient corresponds to a weighted
value in the zone of water table
fluctuations.
-------
Table 13. SATURATED FLOW SUBMODEL
Strengths
Weaknesses
Comments
U)
Ox
Submodel calculates the
steady-state position (and
corresponding flow tube
volumes) of stream lines
extending from the water
table to the drain, using
potential flow theory.
Submodel calculates the mean
travel times for flow
through each tube.
Submodel simulates flow
through a homogeneous, iso-
tropic aquifer.
Two-dimensional flow simu-
lation includes flow to
subsurface, tile drains.
Stream lines occupy fixed
positions as opposed to the
variable situation in the
field.
The travel times are based on
the computed Darcy velocity
divided by the porosity. The
porosity is based on volu-
metric considerations and not
on cross-sectional areas.
Real aquifers are not homo-
geneous and isotropic.
Three-dimensional flow and
radial flow to wells are not
included.
This approximation probably is
suitable in light of other uncer-
tainties in the saturated zone.
The Dupuit-Forchheimer assumptions
are not invoked.
Research in this area has indi-
cated that use of the porosity
yields good results (within a few
percent) when compared to actually
measuring areas of flow. This is
true, even for very tight, clay
soils.
A submodel capable of simulating
a multilayered aquifer would be
desirable at some future time.
The basic techniques are appli-
cable to any geometric flow system.
Analytical solutions are desirable
as they are cheaper to use with
the computer. However, numerical
solutions could be used to simu-
late very complex problems.
-------
Table 14. SATURATED CHEMISTRY SUBMODEL
Strengths
Weaknesses
Comments
Simulates inorganic chem-
istry and denitrification
in the saturated zone.
Flow through each stream
tube is considered to be
via complete piston
displacement.
Mechanical dispersion is
simulated by including
multiple-tube elements.
Inorganic chemical reactions
at each point in space.
Submodel has many options
such as salt injection at
the barrier, simulation of
unsaturated chemistry, etc.
Acid soil chemical reactions
not included.
Denitrification submodel is
very crude.
Flow in the field does not
strictly follow this pattern.
Mechanical dispersion is not
not modeled directly.
Contact times near the drain
may be too small to allow
chemical equilibrium to be
achieved.
User may not be completely
aware of all the options.
These deficiencies will be cor-
rected in the future.
Verification data indicate this
approximation is suitable.
Verification data indicate model
approximation is suitable.
Other uncertainties probably
overshadow this approximation.
The Saturated Chemistry user's
manual describes the options.
-------
Table 15. DRAIN EFFLUENT SUBMODEL
Strengths
Weaknesses
Comments
Combines water from each
flow tube to form the drain
water quality mix.
Water quality mix is not
reacted chemically and could
be unstable from a chemical
standpoint.
Submodel requires consider-
able computer memory for runs
involving many years of
simulation.
This program could be redone using
a different programing approach
to remove this limitation.
OJ
oo
-------
SECTION VI
SUMMARY OF REQUIRED INPUT DATA
From the preceding model description, it is apparent that many physical
and chemical parameters are required for input to the model. Most of
this information is routinely collected or developed in planning and
designing an irrigation project. All of the required field data are
listed below:
DRAINAGE DESIGN
The geometry of the drainage system must be specified including:
a. Drain spacing.
b. Depth to drains.
c. Depth to barrier.
d. Gravel envelope size or design, used to determine the effective
radius.
SATURATED MATERIAL PROPERTIES
Saturated material properties should correspond to values used in design-
ing the drainage system and consequently represent average values for
the assumed homogeneous system. These properties include:
a. Saturated hydraulic conductivities (permeabilities).
b. Porosity.
c. Specific yield or storage coefficient.
UNSATURATED MATERIAL PROPERTIES
In order to simulate unsaturated flow, the following material proper-
ties must be known:
a. Hydraulic conductivity (permeability) as a function of moisture
content.
b. Moisture release curves (tension head versus moisture contents).
Direct measurement of the unsaturated hydraulic conductivity by field
or laboratory methods is difficult, time consuming, and expensive. A
39
-------
limited number of such determinations have been reported in the liter-
ature. Consequently, in lieu of direct measurements, an approximate
procedure known as Millington-Quirk Method (9) can be used.
This approach recognizes that unsaturated conductivity is related to
pore-size distribution. The moisture release curve is used in conjunc-
tion with a semiempirical equation. To assure that the computed and
measured conductivity agree at least at one point, a matching factor
is used. For practical purposes, it is computed as the ratio of the
measured to the computed conductivity at saturation. The technique is
simple enough to allow the determination of the spatial variability of
conductivity.(15) In using the method, moisture contents should be
determined for tensions between 0.05 and 15 atmospheres. Enough points
should be used to adequately determine the shape of the curve.
SOILS ANALYSIS
Master site soils data should provide the following chemistry data:
Soil extract data:
a. Cations (meq/1 of calcium, magnesium, and sodium).*
b. Anions (meq/1 of carbonate, bicarbonate, chloride, and sulfate).
c. Extract ratio of water to soil.
d. If nitrogen is considered, the meq/1 of ammonium, nitrate, and
urea.
Additional data;
a. pH.
b. Cation exchange capacity in meq/100 g of soil.
c. Gypsum content in meq/100 g of soil.
d. Indication if soil is calcareous (line present) or not.
e. Bulk density in g/cm3 of soil.
f. If nitrogen is of interest, the organic nitrogen content in
micrograms/g of soil and the corresponding carbon-nitrogen ratio;
* In this report meq denotes mi11equivalents, meq/1 denotes milliequiv-
alents/liter and mg/1 or mgm/1 denotes milligrams/liter. Also cm^ is
cubic centimeters, cm is centimeters, and g is grams.
40
-------
for use in the transition-state nitrification model, nitrifier-
population salt response relationship.
Normally, these data are collected in 12-inch intervals or some other
suitable scheme down to the barrier. The unsaturated chemistry pro-
gram accepts these data horizon by horizon while the saturated chem-
istry program uses a single or average analysis.
CROP INFORMATION
The different types of crops grown in the area to be studied should be
determined and similar crops grouped. If predictions are to be made
for the entire area, the associated acreage for each group should be
available. For each group, the following information must be obtained:
a. Growing season.
b. Root zone depth and root distribution within the zone.
c. Evapotranspiration rates (if irrigation scheduling is used,
solar radiation and temperature; planting, cover and harvest dates;
cutting dates if alfalfa is included; cutoff date beyond which irri-
gations are not allowed; thaw and freeze dates for the soil; initial
soil moisture in root zone on the thaw date; and the percent soil
moisture in the root zone that may be used before an irrigation is
required).
d. Plant uptake of nitrogen (if nitrogen is considered).
WATER APPLICATIONS
The following data are needed to characterize the water volumes enter-
ing the soil-aquifer system:
a. Irrigation schedules and amounts (for model purposes only that
component entering the soil system is used).
b. Days on which precipitation occurs and the corresponding amounts
(for model purposes only that component entering the soil system is
used).
c. Effective precipitation relationship.
FERTILIZER INFORMATION
Additions of chemicals by fertilizer applications may be included. If
they are, necessary inputs include:
41
-------
a. Fertilization schedules.
b. Amounts of constituents in chemical fertilizers (amounts in
Ib/acre of nitrate, ammonium, urea, calcium, sulfate, and
carbonate).*
c. Depth if plowed in or surface application indicated for each
application of a chemical fertilizer.
d. For each organic nitrogen application the depth of plow down
(no surface application allowed), carbon-nitrogen ratio, and the
application rate in Ib/acre.
IRRIGATION WATER ANALYSIS
The chemical analysis of the irrigation water in meq/1 of NHi4+, N03~,
Ca++, Na+, Mg++, HC03~, Cl", C03=, and S
-------
SECTION VII
INTERFACING WITH CONJUNCTIVE USE MODEL
The return flow quality simulation model described in Volume V is
designed to provide, among other outputs, water flow and quality
predictions in a tile drain or stream. Combination (mixing) of
two or more drains or streams is currently beyond the scope of the
model. This capability could easily be added in the future. The
conjunctive use model (Volume III) can supply this additional mixing
or nodal capability provided a suitable model interface is done.
The soil (unsaturated) portion of the conjunctive use model can be
bypassed in favor of the more sophisticated unsaturated flow and
chemistry programs available in the Volume V model. Also, if the
user wishes to simulate nitrogen transformations, he must utilize
Volume V.
The "tank" type aquifer simulation in the conjunctive use model can
be bypassed in a similar manner in favor of the more detailed satu-
rated flow and chemistry simulations available in the Volume V model.
The user must decide which model or model combinations best fit his
particular application. Applications exist in which one entire model
should be applied, while others would best be handled by some com-
bination of models.
The details of interfacing these two models are left to the user.
The problem is basically one of converting output files to the current
format and time frame for input to the other model. Small interfacing
programs such as the Interface for Chemistry Program will be needed.
The user should consult the appropriate user's manuals for details
concerning the proper input/output files.
43
-------
SECTION VIII
USER'S MANUAL - OVERALL MODEL
INTRODUCTION
This manual discusses details and procedures common to the model
as a whole. Data interfacing between subprograms, bypass options,
and general tips on model application and operation are given.
The model is fully operational on the Bureau of Reclamation CDC
CYBER 74-28 computer located in Denver, Colorado. This computer
is highly time-share oriented with permanent disk storage avail-
able to the user. On this type system, linkage between subpro-
grams within the model can be done easily and the entire model
operated from a remote terminal. The model can also be run on
batch computer systems, but more user effort is required.
SUBPROGRAM LINKAGE
After the basic input data have been entered into the computer,
interfacing between subprograms' is accomplished via binary data
files. These files generally are stored either on magnetic tape
or on permanent disk storage, if available. Output and input for-
mats have been structured to allow direct interfacing without user
intervention. The user need only direct the particular submodel
to read the appropriate binary file(s) created by the previous pro-
gram. Specific discussions of each file are included in the appro-
priate user's manual. Should the user find it necessary to examine
these binary files, a suitable program can be written to read and
print them.
BYPASS OPTIONS
A user may not have the data or require the detail associated with
some models subprograms. The following discussion summarizes bypass
options which have shown utility.
Irrigation Scheduling
The irrigation scheduling submodel can be bypassed if the user
inputs pertinent data to the Unsaturated Flow and Drainout Pro-
grams. This includes crop evapotranspiration data (p. 79) and irri-
gation and precipitation dates and amounts in the case of the Unsat-
urated Flow Program (p. 77). An option (CR0P, p. 76) in the Unsatu-
rated Flow Program allows estimation of crop consumptive use by the
Blaney-Criddle method but is limited to a few crops programed into
the model. The Drainout Program requires data on deep percolation
44
-------
amounts from precipitation and irrigation. These can be supplied
from the Irrigation Scheduling Program or some other source.
Unsaturated Flow And Unsaturated Chemistry
These two submodels tend to be large consumers of computer execu-
tion time. However, they are highly sophisticated programs which
consider many details (including nitrogen transformation and inor-
ganic chemical under nonsteady- and steady-state flow conditions).
The unsaturated chemistry submodel requires input from the Unsatu-
rated Flow Program or a similar source if nonsteady flow is to be
considered. An option in Unsaturated Chemistry (ITEST=1, p. 143)
allows input of steady-state flow and moisture content data from
cards. In this case, additional flow and soil moisture content
data are not needed.
Bypassing of both the unsaturated flow and unsaturated chemistry
submodels can be done quite easily via card input to the Saturated
Chemistry Program. In this case, both quantity and quality of
waters must be specified. The user either enters values for the
volume and quality of the deep percolation water or specifies sim-
ilar values for water entering at the soil surface together with
the onfarm irrigation efficiency. In the latter case, the pro-
gram adjusts both quantity and quality in proportion to the
efficiency.
There is no direct option which allows input of data to the Satu-
rated Chemistry Program from the Unsaturated Flow Program. If
Unsaturated Flow is utilized, Unsaturated Chemistry must be used
to generate input to Saturated Chemistry. Another program could
be constructed to take unsaturated flow output from Unsaturated
Flow, assign qualities to these flows, and punch a card deck suit-
able for input to the Saturated Chemistry Program.
The user should note that a bypass of Unsaturated Chemistry bypasses
the nitrogen transformations in the unsaturated zone. There is no
alternate means of simulating these reactions within the model. It
is possible to use the Saturated Chemistry Program to simulate salt
reactions in the unsaturated zone. In this case, some characteristic
leaching moisture content(s) is specified in place of the saturated
value(s). Each tube in Saturated Chemistry then represents a spe-
cific soil profile and each tube element becomes an unsaturated soil
segment.
Additional options in Unsaturated Chemistry allow the nitrogen trans-
formations and salt reaction portions to be individually bypassed.
45
-------
Thus, it is possible to consider nitrogen transformations and move-
ment without the salt reactions. Conversely, salt reactions and
movement can be simulated without the nitrogen transformations.
Saturated Chemistry
This submodel is not needed if the simulation is terminated at the
water table and the program is not being used to simulate chemistry
in the unsaturated zone. It is possible to use Saturated Chemistry
as purely a routing-type model with no salt chemical reactions con-
sidered. In this case, SBYPAS (p. 183) is set equal to one. This
option is extremely useful in cases where only nitrogen species are
being considered in the saturated zone. If IDN (p. 183) is set equal
to one, a denitrification subroutine is called from Saturated Chemistry.
If the parameter equals zero, the subroutine is bypassed.
Another option in Saturated Chemistry allows injection of salts at
the barrier. In this case IINJE (p. 183) is set equal to one, and addi-
tional data (Card Group 4, p. 186) must be inserted into the input deck.
This option allows salts to be added into the system from a source
such as a shale bed. Pick's Law is utilized.
Still other options in Saturated Chemistry allow user inputed values
for ion exchange coefficients and CC>2 partial pressures. Appropriate
sections of the user's manual should be consulted for details.
Drain Effluent Prediction
The submodel can be bypassed any time there is no need to mix waters
from the various saturated flow tubes. This can occur when the
simulation is terminated at the water table, when drain discharges
without any associated qualities are being simulated, or when effects
on the aquifer are desired independently of drain quality.
The Drain Effluent Prediction Program must be utilized if quality
predictions at the drain are desired.
Saturated Flow
The Saturated Flow Program can be bypassed if values for flow tube
volumes and top widths can be determined from alternate sources.
These could be estimates or other simulation models. The program
is not needed if only monthly drain discharges with no associated
qualities are the objective. In the latter case, these values are
output from the Drainout Program.
46
-------
TIPS ON MODEL APPLICATION
Use Of Bypass Options
Selection of the total model or one or more bypass options depends
on many factors. The user must determine his needs based on data
resources, time and funding allocations, desired accuracy, and
model sensitivity as applied to a particular situation. No fixed
rules have yet been developed, but common sense applied together
with a critical evaluation of the problem should help provide the
necessary guidance.
The flexibility built into the model was placed there to allow a
wide range of applications. However, care should be taken not to
abuse this feature. A program sequence used in one application may
be entirely inappropriate for another.
Several model configurations have shown utility based on experience
over the past few years. These are illustrated in Figure 2 and in
Figures 6 through 11. Additional sequences are possible based on
the model bypass options available to the user. The sequences shown
illustrate some of the more common applications possible with this
model. Figure 2 contains all of the subprograms contained in the
model. This sequence allows the maximum amount of detail which the
model can supply. This configuration is recommended whenever the
data and study requirements indicate that a maximal effort be put
forth. It is not recommended for screening or reconnaissance-type
applications.
The configuration shown in Figure 6 is a useful computer time-
saving sequence when the unsaturated zone can be ignored. This
can occur when relatively "clean" surface soils are encountered,
nutrients are not being considered, and "ball park" type answers
are wanted for large acreages. With this sequence, many computer
runs can be made at a low cost. No information can be obtained
concerning chemical changes in the unsaturated zone (e.g., the
location of salt precipitation within the unsaturated zone cannot
be determined).
Figure 7 contains a configuration which is similar to that in Fig-
ure 6 except that piston flow at a constant 'leaching1 moisture
content is considered in the unsaturated zone. The Saturated Chem-
istry Program is used to simulate salt reactions (no nutrient trans-
formations) in the unsaturated region. It generates output which
can again be inputed to the same program but at the saturated mois-
ture content. This arrangement uses only slightly more computer
time than that in Figure 6 and includes simulation in the unsat-
urated zone.
47
-------
IRRIGATION
SCHEDULING
DEEP
PER. VOLUMES
r
I
DRAINOUT
DEEP
PER. VOLUMES
SATURATED
FLOW
DRAIN
EFFLUENT
PREDICTION
WATER APP /
DEEP PER.VOL.
& QUALITY
SATURATED
CHEMISTRY
Assumes chemistry and flow in the unsaturated zone have
insignificant effects on quality of drain effluent. Irrigation
scheduling may be utilized or by passed.
FIGURE 6. MODEL BY PASS CONFIGURATION NO. I
48
-------
IRRIGATION
SCHEDULING
DRAINOUT
WATER APR
&
QUALITY DATA
DEEP
PER. VOLUMES
SATURATED
FLOW
DRAIN
EFFLUENT
PREDICTION
SATURATED
CHEMISTRY
SATURATED
CHEMISTRY
More crudely approxmates unsaturated flow by assuming leaching
fakes place by piston displacement at some constant leaching
moisture content.
FIGURE 7. MODEL BY PASS CONFIGURATION NO. 2
49
-------
MONTHLY
DRAIN
DISCHARGE
WATER APR
&
QUALITY DATA
TUBE VOLUMES
&
TOP WIDTHS
TUBE TRAVEL
TIMES
SATURATED
CHEMISTRY
SATURATED
CHEMISTRY
DRAIN
EFFLUENT
PREDICTION
All water flow parameters must be estimated or obtained from
another source. Chemistry routines utilized whenever possible.
FIGURE 8. MODEL BY PASS CONFIGURATION NO. 3
50
-------
IRRIGATION
SCHEDULING
TUBE VOLUMES
TOP WIDTHS
TUBE TRAVEL
TIMES
MONTHLY
DRAIN
DISCHARGES
UNSATURATED
FLOW
INTERFACE
UNSATURATED
CHEMISTRY
SATURATED
CHEMISTRY
DRAIN
-EFFLUENT
PREDICTION
Considers salt and nitrogen chemistry in unsaturated zone. Groundwater
flow parameters estimated
FIGURE 9. MODEL BY PASS CONFIGURATION NO. 4
51
-------
IRRIGATION
SCHEDULING
UNSATURATED
FLOW
INTERFACE
UNSATURATED
CHEMISTRY
Detailed unsaturated zone
FIGURE 10. MODEL BY PASS CONFIGURATION NO. 5
-------
IRRIGATION
SCHEDULING
UNSATURATED
FLOW
DRAINOUT
Allows design of drainage system using either output from
Unsaturated Flow or deep percolation values from
Irrigation Scheduling.
FIGURE II. MODEL BY PASS CONFIGURATION NO. 6
53
-------
Figure 8 illustrates a sequence which is useful when flow through
large aquifers to open drainage channels such as stream and rivers
is being modeled. In this case, the user inputs estimates for all
the flow parameters and utilizes a simplified chemistry sequence.
Nitrogen transformations in the unsaturated zone cannot be modeled.
The configuration in Figure 9 is somewhat similar to the one in
Figure 8 but allows nutrients as well as salts to be considered
when groundwater flow parameters must be estimated. In this
case the unsaturated flow subprogram must be included for nonsteady
flows, and the irrigation scheduling subprogram is recommended.
Figure 10 illustrates an abbreviated sequence which allows maximum
detail in the unsaturated zone. This is useful when the effects
over time of various cultural practices are desired. Effects of
high irrigation efficiencies, gypsum applications, high salinity
and/or high sodium irrigation waters, etc., can be tested in some
detail in the root zone.
The abbreviated configuration shown in Figure 11 allows design of
drain spacings given porous media properties, depth to barrier,
effective drain radius, depth to drains, and the deep percolation
schedule. Irrigation scheduling can be used to input monthly deep
percolation volumes into drainout, thus bypassing unsaturated flow.
Intermediate Output
The model is composed of several subprograms which have been inter-
faced. Thus, intermediate results are generated at various stages
during a simulation run. The user should take particular care to
examine some of these results (not necessarily the entire binary
interface files) and make certain he is satisfied before proceed-
ing to the next step. Often an error has been made somewhere in
the input data or assumptions. If an error is present, early detec-
tion could avoid a costly rerun of some or all of the model. Check
the results (intermediate and final) against your experience and
judgment and any verification data which may be available. Never
trust results merely because the computer has generated them. Gar-
bage (data and assumptions) into a machine always equals garbage
out, and good inputs and assumptions yield qualified predictions.
Printed Output Volume
The output options in the model allow the generation of large vol-
umes of printed output. Although a certain amount is needed to
verify intermediate results and list the final answers, an attempt
should be made to minimize the volume. Output can be dumped to
magnetic tape or microfiche in lieu of printing, but again the
question should be asked, "Is all this really necessary?"
54
-------
Accuracy Versus Computer Run Times
Within the model, several parameters 'are specified which allow the
user to increase the accuracy of the results at the expense of com-
puter execution times. In most cases, a range of values is listed
to give the user an idea of what numbers to select. The user must
use his own judgment as to what accuracies are demanded in his par-
ticular application. Some experimentation may be necessary to more
closely refine some of these criteria.
Conversion To Other Computer Systems
The computer model as listed in Volume V is operational on a CDC
CYBER 74-28 machine under the FORTRAN Extended Compiler. Conver-
sion to other CDC 3000, 6000, and 7000 series computers should be
possible with a minimum of problems. Conversion to machines such
as IBM and Univac is possible, but a number of difficult conver-
sion problems can be expected. The user should use the test run
included in this volume as a check on the validity of his partic-
ular conversion.
Some consideration should be given to utilizing this model on the
Bureau of Reclamation computer system. A remote terminal can be
tied into the system and used to run the model with a minimum of
effort. In addition, Bureau personnel can readily assist with any
program difficulties.
55
-------
SECTION IX
USER'S MANUAL
IRRIGATION SCHEDULING PROGRAM
INTRODUCTION
This program computes irrigation schedules using daily historic pre-
cipitation, temperature and solar radiation data, and associated crop
and soil data. Crop evapotranspiration is computed using the Jensen-
Haise procedure (7,8). Irrigation dates and amounts are determined
on the basis of a soil moisture budget accounting procedure.
If desired, a file containing crop evapotranspiration data and pre-
cipitation and irrigation applications can be written in a format
suitable for input to the Unsaturated Flow Program via an interface
utility routine.
Likewise, a file containing deep percolation amounts from rainfall
and irrigation can be written in a format suitable for input to the
Drainout Program via an interface utility routine.
For comparative purposes, weighted-average, monthly irrigation infor-
mation can be computed from the irrigation events for each crop. These
values can then be compared with a monthly summary of irrigation events
computed using the weighted evapotranspiration values from all crops.
The program reads on file 2 (input) and writes on files 3, 4, 5, and 6
(output). The FORTRAN program requires about 240,000 octal words on
the Bureau computer. The storage locations in the computer must be set
to zero before loading this program.
DATA INPUT CARD SEQUENCE AND DATA LIMITATIONS
CARD GROUP 1
Column Var Var Variable
No. type* name description
Card 1 of 5
1-5 A IZIP Zip code for study area.
6 X - Region number.
Note: A is alphanumeric.
R is real (floating point)
I is integer (fixed point)
X is blank or space.
56
-------
Column Var
No . type
7-36
37-42
45
(The following
46-48
49-50
51-53
54-55
79-80
Card 2 of 5
1- 6
7-72
79-80
Card 3 of 5
1- 6
14-15
A
R
I
10
I
I
I
I
X
X
A
X
X
I
Var
name
ILOCA
ELEV
NELFLG
columns not
LATDEG
LATMIN
LONGDG
LONGMN
-
-
IUSBR
-
-
NOYRS
Name an<
Elevati<
Flag foi
T AT>
-------
Column
No.
30
31-35
36-40
52
54
56
Var Var Variable
type name description
I NCPTS Number of combinations of weighted crop
percents.
Min = 0 (if NCPTS = 0, do not include
Card 4).
Max = 5.
R PCTWPB Percent of winter snow that enters the
R PCTWPC soil at the time of soil thaw. The
remaining percent of the snow is con-
sidered runoff.
Note: Winter precipitation has daily
evaporative demand subtracted; the
remaining precipitation is summed
throughout the winter (carried over
from 1 year to the following year)
and the accumulated value used at
soil thaw.
PCTWPB is percent for bare ground.
PCTWPC is percent for ground with cover
such as alfalfa or pasture.
I IPUNCH Flag for writing precip, irr, and ET on
file 4 in a format suitable for input
to the unsaturated flow program.
(Interface utility.)
If card output is desired, cards may be
punched from file 4.
0 = No write.
1 = Write.
I IFLGCL Flag for printing climatic data (file 3).
0 = No print.
1 = Print.
I IFLGIR Flag for printing irrigation schedules
by year and by crop.
0 = No print.
1 a Print.
Note: If both IFLGCL and IFLGIR are zero,
only a summary will be printed which
shows annual values for all years for
each crop.
58
-------
Column
No.
58
Var Var
type name
I IWIDTH
60
ICARYO
62
IWTDP
64
IWTIR
66
IWTIRD
Card 4 of 5
79-80
1- 6
X
X
Variable
description
Flag for width of printout of irrigation
schedules.
0 = 132-column printout.
1 » 72-column printout.
Flag for carrying over soil moisture
and winter precipitation to the follow-
ing year.
0 = No carry-over.
1 - Carry-over.
Note: If ICARYO = 1, NDB and NDE should
be 1 and 366, respectively.
Flag for writing deep perc from rain and
irrigation in format for input to
Drainout Model.
0 = No write.
1 = Write.
Flag for computing weighted monthly
irrigation parameters.
Note: This weighted average computed by
irr event was done to compare weighted
et values from all crops. An annual
summary is written on file 6 when this
flag is set.
0 = Annual summary not printed.
I = Annual summary printed.
Flag for printing detailed event by event
monthly summary of irrigation weighted
information on file 6.
0 » Detailed summary not printed.
I = Detailed summary printed.
Card number.
Zip and region.
59
-------
Column
No.
16-45
Var
type
R
Var
name
CPCT
Variable
description
The percentage of each crop acreage to the
total acreage, as it appears chronologi-
cally in card 9, three columns per value,
maximum of 10 values (10 crops) .
Card 4 is repeated for each combination of weighted crop percents.
There should be the same number of these cards as NCPTS.
Card 5 of 5
1-15 X -- Zip, region, and other ID.
16-40 I IPAR Flag for type of irrigation parameter to
be weighted and printed in monthly for-
mat. If IPAR = 0, do not do it. If
IPAR = 1, do it. Place the flag in the
proper column, see following list.
Col parameter:
16 Evapotranspiration
18 Precipitation
20 Precip runoff
22 Precip entering soil
24 Precip deep percolation
26 Soil moisture added by precip
28 Net et (et reqt to be met by irr)
30 Total irrigation application
32 Drift loss from sprinkler
34 Runoff loss from gravity irr
36 Irr appl infiltrating the surface
38 Soil moisture added from irrigation
40 Deep percolation from irrigation
Note: Parameters listed in col 16-28
Card 5 is included only if IWTIR = 1.
CARD GROUP 2
Card 1 of 5
1-6 X - Zip and region.
9-10 I NTHB Number of the beginning month of poten-
tial evapotranspiration calculations.
60
-------
Column
No.
10-14
Var Var
type name
I KPCF
15-19
20-21
22-23
I IDATA
I ICYCLE
24-25
NCROP
26-29
30-33
NDETB
NDETE
Variable
description
Flag for computing et using crop curves
from cards or polynomials in the CPCOEF
subroutine.
0 = Read crop curves from cards 12 and 13.
1 = Compute crop curves from polynomials.
Columns used for additional indentifi-
cation - farm no, field no, etc.
Flag for computing irrigation schedules.
0 = Do not compute schedules.
1 = Compute schedules.
Flag for reading another set of crop data
(card 9) or reading climatic data (cards
1-8).
0 = Read card 1 and new year of dim data.
1 = Read card 1, NYR and IREAD only.
Note: This prevents the inclusion of the
rest of data in card 1.
2 = Read card 9 for another crop.
Crop identification number:
1 = Small grain
: Beans
•• Peas
= Potatoes
= Sugar beets
= Corn and sorghum
= Alfalfa
= Pasture
Weighted crop does not use a crop
curve but requires a crop number to
operate properly.
Note: For weighted crop use any crop num-
ber (must be .GT. 0),
Julian date (1-366) of planting or
beginning of growing season.
Julian date (1-366) of harvest or end
of growing season.
61
-------
Column
No.
34-37
Var Var
type name
I NDCOV
38-41
43-60
I NCUTNG
A ICROP1
Variable
description
Julian date (1-366) of effective crop
cover date.
Note: The Julian dates used should be
for a leap year (366).
The program uses a 366-day array but
leaves out Feb 29 if the year is not
a leap year. Dates for southern Idaho
conditions are listed as follows (after
Jensen, Wright, and Pratt; estimating
soil moisture from climate, crop, and
soil data; Transactions ASAE, Vol 14,
No. 5, Sept-Oct, 1971, pp. 954-959,
unpublished Appendix A, May 1971):
Small grain At heading
Beans Bloom or about
50 days after
planting.
Full bloom or 70 days
after planting.
About 80 days after
planting.
About 110 days
after planting.
About 10 days after
tasseling on corn
and heading for
sorghum.
All season except
30 days after
growth begins in
the spring and
20 days after
cuttings.
Note: Cover data for alfalfa are not
required, as they are taken care of
internally in the program.
Pasture All season except
30 days after
growth begins in
the spring.
Number of cuttings for alfalfa.
Eighteen spaces available for crop name.
Peas
Potatoes
Sugar beets
Corn or sorghum
Alfalfa
62
-------
Column Var Var
No . type name
61-66 A IPIANT
67-72 A ICOVER
73-78 A IHARV
79-80 X
Card 2 of 5
1-8 X -
9-12 I NCUT
13-18 A ICIJT
Variable
description
Plant data in month/day nomenclature.
Cover date in month/day nomenclature.
Harvest date in month/day nomenclature.
Card number.
Zip, region, and year.
Julian date (1-366) of first cutting
for alfalfa. Use this card only for
alfalfa.
First cutting date in month/day terms.
Subsequent NCUTS and ICUTS are entered on this card as 14, A6 up
to a maximum of seven cuttings which ends in column 78.
79-80 X
Card 3 of 5
1-20 X
21-75 R RP
Card 4 of 5
1-25 X
26-75 R RD
Card number.
Crop type and other identification.
Crop coefficients from the Jensen-Haise
crop curves at 10 percent intervals,
starting at 0 and ending on 100 percent
cover (a total of 11 values, 5 columns
per coefficient) .
Crop type and other identification.
Crop coefficients at 10-day intervals
beginning with the 10th day after full
cover and ending with 100th day (a
total of 10 values, 5 columns per
coefficient).
Note: Include cards 3 and 4 only if KPCF = 0 on card 1.
63
-------
Column
No.
Card 5 of 5
1-15
16-20
21-25
26-30
Var Var
type name
X
I
IBGSEQ
IENSEQ
JENDAY
31-35
36-40
R
AVAIL
DSMBG
41-45
R
ADEPCT
46-50
R
ADEPC1
51-55
FMEFF
Variable
description
Zip, region, and other ID if desired.
Julian date (1-366) of soil thaw. Soil
moisture budget computations are begun on
this date.
Julian date (1-366) of soil freeze. Soil
moisture budget computations are stopped
on this date.
Julian date (1-366) of last allowable
date of irrigation in the season when
properly specified allows a drying period
prior to fall harvest.
Inches of available soil moisture in the
root zone of the crop specified in
card 1.
Depleted soil moisture of soil moisture
budget computations.
If ICARYO = 1 on card 1C only the DSMBG
value for the first year will be used.
Allowable depletion percentage. This is
the percent of the available soil
moisture which the crop is allowed to
use before an irrigation is called for.
Allowable depletion percentage for the
first irrigation. Oftentimes this
percent is less than ADEPCT because
the root zone is not fully developed
by the first irrigation.
Irrigation efficiency (in percent) for
the type of system involved. Irrigation
efficiency is defined as soil moisture
added to the root zone divided by water
delivered, expressed as a percent.
64
-------
Column Var Var Variable
No- .type name description
56-60 R DRFTLS Drift loss (in percent) if a sprinkler
system is involved. This percent would
include evaporation loss and is a per-
cent of the water delivered to the field,
not the water applied to the surface.
61-65 R ROFLS Runoff loss, in percent, of the water
delivered to the field.
66-68 X
69-70 A IRRT Irrigation type.
GR - Gravity
Sp - Sprinkler
Repeat cards 1-5 for each crop. Remember, that on the last 1 card
in group 2, ICYCLE is set to 0 to begin reading card group 3 and
climatic data.
CARD GROUP 3
Card 1 of 32
1-6 X - Zip and region.
7-8 I NYR Last two digits of the year.
To terminate the job set NYR = 0 by
using a blank card.
13-14 I IREAD Flag for reading climate data.
0,1 = Read only precipitation, temper-
ature, and solar radiation, in
that order.
2 = Read precipitation, evaporation,
temperature, and solar radiation
in that order.
3 = Read precipitation, relative humid-
ity, evaporation, temperature, and
solar radiation in that order.
4 = Read precipitation, wind, relative
humidity, evaporation, temperature,
and solar radiation in that order.
5 = Read only ETP (potential evapotran-
spiration) that has been computed
previously and placed on cards.
65
-------
Column
No.
Var Var
type name
15-18
19-22
NDB
NDE
23-30
ATERM
31-38
39-47
48-59
R
Variable
description
A minus sign - in front of any of the
above IREAD numbers prevents the program
from using any crop information.
A minus sign is required for the first
year of climatic data. For all the
following years no minus is used.
Julian date of beginning climatic data.
Julian date of ending climatic data. If
a full year data are processed, NDE
should be 366. The program operates on
a 366-day array and leaves out Feb 29 if
the year is not a leap year.
Average daily maximum temperature for
the warmest month or the CT coefficient
in the Jensen-Haise equation,
ETP = CT*(T-TX)*RS.
If the temperature is input, it must be
preceded with a minus sign as a flag.
The minus causes CT to be computed in
subroutine JHCOEF.
Average daily minimum temperature for the
warmest month or the TX coefficient in
the Jensen-Haise equation. A minus sign
is not needed with this temperature
input. See climatic summary of United
States for ave daily max and min temps.
Conversion from gram-calories per square
• centimeter (Langleys) to inches of
evaporation equivalent. CONVRT is used
to convert solar radiation in Langleys
to inches.
CONVRT = 0.000673 if solar radiation
is in Langleys.
CONVRT = 1.0 if solar radiation is
in inches.
SUMETP, These variables are read in as zeros to
SUMEVP start sums.
BTERM
CONVRT
66
-------
Column
No.
73-78
79-80
Var Var
type name
Z OUT
Variable
description
This value read in as a blank for
printing blanks in climatic data
output.
Card number.
The following description applies to cards 2 through 32 which are
climatic data - set up for 1 year on 31 cards - 1 card for each day
of the month and 12 columns each representing a month for that
particular day.
Card 2 of 32
1- 6
7- 8
9-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
X
X
I
X
R
R
R
R
R
R
NDAY
DATA
DATA
DATA
DATA
DATA
DATA
Zip code and region.
Last two digits of the year.
Day of the month (1-31).
Identification for data type:
Datatype
(Col 11-15)
Prep
Wind
Rel H or
Dewpoint
Evap
Temp
RS (Solar)
ETP
Card No.
(Col 79-80)
2
3
4
5
6
7
8
Data for January corresponding to NDAY.
Data for February.
Data for March.
Data for April.
Data for May.
Data for June.
67
-------
Column Var Var Variable
No. type name description
46-50 R DATA Data for July.
51-55 R DATA Data for August.
56-60 R DATA Data for September.
61-65 R DATA Data for October.
66-70 R DATA Data for November.
71-75 R DATA Data for December.
79-80 R DATA Card number.
Detailed information on data input is as follows:
Prep Precipitation in inches for the day - input as
hundredths of an inch.
Wind Wind in miles per day - no decimals.
Rel H Relative humidity in percent - no decimals.
or
Dewpoint Dewpoint in degrees F.
Evap. Pan evaporation in hundredths of an inch
(comparative purposes).
Temp Temperature input as maximum daily temperature in
the first three columns of the five-column field,
and as minimum daily temperature in the last two
columns of the five-column field. For example:
Field XXXXX
Temp 9852
Maximum daily temp = 98
Minimum daily temp = 52
Average daily temperature may also be used. In
the case where negative max or min temps are
involved, the average temp must be used.
(Note: A negative average may be used.)
68
-------
Column Var Var Variable
No. type name description
RS Solar radiation in Langleys or inches - see CONVRT
on card 1. Langleys - no decimals, inches -
hundredths.
ETP Potential evapotranspiration computed externally
and input for scheduling purposes.
The last data card is a blank card to exit the program - See NYR
on card 1.
69
-------
SECTION X
USER'S MANUAL - UNSATURATED FLOW PROGRAM
INTRODUCTION
This program is designed to simulate one-dimensional unsaturated flow
between the ground surface and the water table. It includes the
infiltration, redistribution, plant root extraction, and drainage of
soil water under a growing crop. In particular, this program was
developed as part of a water quality simulation model with the water
contents and movement within the unsaturated zone as well as the deep
percolation to the water table acting as input to other programs. A
generalized diagram of the system modeled appears in Figure 12.
The program was originally written at the University of Arizona under
a Bureau of Reclamation contract. The theoretical basis has been
described in Dutt, et al. (5). However, a number of program improve-
ments and modifications have been made including:
a. Improved restart ability, including use of tape, cards and
printed output.
b. Addition of ability to simulate more than 1 year with a single
run, including commencement of computations with a thaw date and
stopping them on a freeze date.
c. Improved method of handling the upper boundary condition when
there is no water at the surface and the zero flow condition
prevails.
d. Addition of an option to compute the thaw date application when
the profile is assumed at field capacity and multiple years are run.
e. Improved computation of the consumptive use to assure mass bal-
ance is maintained.
f. Consolidation of all references to a specific soil in program
M0ISTRE by calling a new subroutine named PR0P.
g. Development of a special form of PR0P in which a particular
functional relationship is used. This permits soils properties for
conductivity and diffusivity to be defined by card inputs, eliminat-
ing the need to modify the source deck.
h. Elimination of certain card groups when their options are not
used.
-------
GROWING PLANT OR CROP
SOIL SURFACE
,J
§
2
S>
©
h4
©
©
NJ
E
©
tn
c
o
O
sz
2
13
>4
•
UJ
Q
Q-2
Q-l
-Top boundary set at
TBC when water at
surface except at start
of run if westart or
group 4 cards read
hINTERNAL NODES
Q = number of nodes
n H0RIZON(0) ,
u " DELX '
(integer arithmetic)
NODE NUMBERS
-BOTTOM BOUNDARY
(stationary, may or may not
represent water table)
-Bottom boundary set at
TBC except at start of run
if restart data or group
4 read
FIGURE 12. MODELED SYSTEM - UNSATURATED FLOW
71
-------
i. Modification of input formats to pack more data on a single
card.
This manual updates and replaces the user's manual found in Dutt,
et al., (5), pp. 31-35.
UNIT NUMBERS
Five input-output devices or files are referenced by the program:
Unit
designation
INPUT
0UTPUT
PUNCH
10
Unit description
or purpose
Card input is accomplished by READ statements
without a unit number. The common input tape
or unit (i.e., the computer system's standard
input device) is understood.
Printed output is generated by PRINT state-
ments without a unit number. The common out-
put tape or unit (i.e., the computer system's
standard output device) is understood.
Records are written to the punched card
device of the computer by PUNCH statements.
The unit to which computed results are written
on a saved magnetic tape. An entirely new
tape may be generated or results can be added
to a previously created tape. This tape is
used as input to the interface program.
The unit to which restart data is written on
a saved magnetic tape.
It is noted that tape 4 (5, p. 32) is no longer used.
DATA DECK STRUCTURE
There are eight basic groups of input cards in the data deck, as shown
in Figure 13, and the variables are defined below.
Groups are numbered to conform with those of the original University
of Arizona user's manual.(1) The consumptive use routine (C0NUSE)
described in that publication does not correspond to the one in our
program. An earlier version, in which the Blaney-Criddle constants
and climatic data for two crops (barley and milo) are built into the
program, is used. Group 5 inputs are not accepted.
72
-------
END OF INPUT DECK
Repeot for eoch year of run
GROUP I
GROUP 2
/TITLE CARD
I FOR APPLIC.
GROUPS
1 ) Denotes card type within a group
START OF INPUT DECK
FIGURE 13. DATA DECK STRUCTURE UNSATURATED FlDW SUBMODEL
73
-------
Group 8 has been added to allow input of soil properties when Campbell's
method is employed to define diffusivity and conductivity relations.
The restart deck does not have a group number, and is described sepa-
rately in a later section.
CARD GROUP 1 - Input-Output Options and Run Conditions
Column
No.
Card 1 of 2
1-5
6-10
11-15
16-20
Var Var Variable
type name description
IPUNCH Option to allow restart data to be
saved at the end of a run. If
IPUNCH = 1, restart data are to
be punched
IPUNCH = 2, restart data are to
be written on tape 10
IPUNCH = any other value, restart
data are not saved
IRESTR Option to allow restart data to be
used as input for this run. If
IRESTR = 1, restart data are read
from cards
IRESTR = 2, restart data are read
from tape 10
IRESTR = any other value, restart
data are not used
ISAVE Flag to specify that a previously
created output tape 5 is to be used.
It is then necessary to position the
tape to the correct position for the
initial write, using INFIL5 and
LLSTRT (below).
If ISAVE = 1, use a previously saved
tape 5.
ITENTH Option to allow printed results to be
outputed at 0.1-day intervals.
If ITENTH = 1, print results every
tenth day.
74
-------
Column
No.
21-25
26-30
31-35
36-40
41-45
Var Var
type name
I INFIL5
LLSTRT
MMST0P
IST0P
IDEF
46-50
R FCAP
51-55
I IP0PT
Card 2 of 2
1-5 I AA
6-10 I BB
11-15
CC
Variable
description
Initial file number on tape 5 for the
first write of output data (see
RESTART OF PROBLEM).
Julian day number of first day for
which computations are made for the
initial year of this run.
Julian day number of last day in the
last year (IST0P) for which compu-
tations are made.
Stopping or ending year for this run.
Option to select the deficit or
depleted moisture method of deter-
mining the thaw date application.
If IDEF = 1, use this option (see
DEPLETED MOISTURE
OPTION).
Volumetric moisture content to be
used as the field capacity value
when the depleted moisture option
(IDEF) is used. FCAP is expressed
as a decimal (dimensionless).
Option to select the special form of
subroutine PR0P to compute conduc-
tivities and diffusivities as a
function of moisture content.
If IP0PT = 1, then use this form (see
UNSATURATED FLOW PROP-
ERTIES) . A Group 8
card must then be
included.
If AA = 1, print the initial input data.
If BB = 1, print computed moisture con-
tents on a daily basis.
If CC = 1, read initial moisture con-
tents from Group 4 cards.
75
-------
Column Var Var
No. type name
16-20
21-25
26-30
I LL
I MM
R BBC
31-35
R TBC
36-40
41-45
I YEAR
I CR0P
Variable
description
Julian day number of starting day in a
full year (i.e., first day for all
years except initial).
Julian day number of starting day in a
full year (i.e., last day for all
years except IST0P).
Volumetric moisture content at the
bottom boundary expressed as a deci-
mal (dimensionless).
If BBC = 1, then the bottom boundary
is set at the TM value.
BBC = 2, then the bottom boundary is
set at the TS value.
BBC = any other value, the bottom
boundary is set at this value.
Volumetric moisture content at the top
boundary expressed as a decimal
(dimensionless).
If TBC = 1, then the top boundary is
set at the TM value.
TBC = 2, then the top boundary is set
at the TS value.
TBC = any other value, the top bound-
ary is set at this value.
The year number corresponding to the
initial year of this run.
Crop selector used in computing con-
sumptive use.
If CR0P = 1, use the Blaney-Criddle
formula with the con-
stants in the program
for barley.
CR0P = 2, use the Blaney-Criddle
formula with the constants
in the program for milo.
CR0P = 3, use consumptive use values
read from cards in Group 7,
76
-------
Column
No.
46-50
56-60
61-65
66-70
71-75
Var Var
type name
I M
R
R
R
DELX
TS
TM
TD
76-80
R
SM
Variable
description
Minimum number of time intervals per
day for applying the finite differ-
ence solution.
Distance between nodes for the finite
difference solution (cm).
Maximum (i.e., usually taken as the
saturation value) moisture content
expressed as a decimal (dimensionless)
Arbitrary moisture content expressed as
a decimal (dimensionless).
Minimum (i.e., usually taken as perma-
nent wilting point) moisture content
expressed as a decimal (dimension-
less). Consumptive use cannot be
withdrawn by the sink term when this
value is reached.
Arbitrary initial moisture content for
internal nodes (i.e., all nodes
except top and bottom boundary) used
when the CC and/or restart option
are not used, expressed as a decimal
(dimensionless).
CARD GROUP 2 - Water Applications
Card 1 of 2
1-5
6-10
11-13
I IC0DE A numeric code that can be used to
identify the crop, climatic condi-
tions, irrigation method, etc. It
is merely displayed on output.
I IYEAR Year number identifying the year of
inputs. It is merely displayed on
output.
I APPS Number of water applications contained
in the type 2 cards of this group
for this year (maximum number is 50),
77
-------
Column
No.
14-77
Var Var
type name
A IDENT
Card 2 of 2
11-12
13-14
15
M0N
I DAT
ADENT
16-20
AMI
Variable
description
Alphanumeric description to identify
the water applications. It is merely
displayed on output and may include
information on the crop, climatic
conditions, and irrigation methods.
Month of water application where
January is 1, February is 2, etc.
Date within month of water application.
An alphanumeric identifier to denote
the source of water, for example.
I = irrigation, R or P = rainfall,
S = snowmelt, etc. It is merely
displayed on output.
Amount of application as a volume per
unit area (inch3/inch2 or inch).
Water is applied at the start of the
day.
CARD GROUP 3 - Soil Identification and Depth
Card 1 of 2
1-2
Card 2 of 2
3-10
11-20
IDENT
Number of soil horizons.
Alphanumeric soil horizon identification
which is merely displayed on output.
R
H0RIZ0N Depth from land surface to lower bound-
ary of horizon (cm).
CARD GROUP 4 - Initial Moisture Contents
Card 1 of N
1-10
11-20
etc.
R T0(I)
Initial volumetric soil moisture content
of each depth node I including top
and bottom boundary nodes expressed
as a decimal (dimensionless). N is
the number of nodes.
78
-------
Column Var Var Variable
No- type name description
CARD GROUP 6 - Root Distribution
Card 1 of 1
1'10 R KP(I) Decimal fraction of roots located in
11-20 the Ith foot of the soil profile
etc. (dimensionless). KP(1) is for the
top foot, KP(2) is for the second
foot, ... KP(6) is for the sixth
foot.
CARD GROUP 7 - Consumptive Use Values
Card 1 of 1
1-5 I IC0DE A numeric code that can be used to
identify the crop, climatic condi-
tions, irrigation method, etc. It
is merely displayed on output.
6-10 I IYEAR Year number identifying the year of
inputs. It is merely displayed on
output.
11-15 R U(I) Semimonthly consumptive use in the
16-20 entire root zone for period I (acre-
etc. inches/acre or inches). 1=1
denotes January 1 through January 15,
1=2 denotes January 16 through
January 31, etc. Thus, the first
period in a month is from the. 1st-
15th, while the second is from the
16th to the end of the month.
CARD GROUP 8 - Unsaturated Flow Properties
Card 1 of 1
1_10 R KSAT Saturated hydraulic conductivity
(cm/day).
H_20 R BEMP Empirical constant b.
21-30 R AIRENT Air-entry potential (cm).
79
-------
NODE SPACING
An idealization of the unsaturated flow system is shown in Figure 12.
The horizon depths given by the group 3 cards, root zone distributions
given by the group 6 card, and the nodal- system are indicated. The
number of nodes, Q, is based on the depth to the bottom of the deepest
horizon, H0RIZ0N (0), and the DELX value of card 2, group 1.
0 - H0RIZ0N(0) ,
V " DELX
The right hand expression is truncated because Q is integer. For
example, if H0RIZ0N(0) = 146. cm and DELX = 5.00 cm
a 29.2+1 = 30.2 -»• 30 [2]
O » \J \J
If H0RIZ0N(0) = 150. cm and DELX = 5.00 cm it follows that
= 31. -> 31 [3]
However, the internal representation of real numbers that are con-
sidered exact on input generally results in a stored value that is
slightly different. Thus, in the previous example, the depth might
be represented internally by a number slightly less than 150.0 cm.
The number of nodes would then be computed as 30 which is a totally
unexpected result.
In general, the use of any real numbers in an integer computation
can produce unexpected results. Since it is necessary to know Q when
initial conditions are read, it is recommended that the H0RIZ0N(0)
depth be adjusted to assure the desired number of nodes. In the pre-
vious example, this is accomplished by setting H0RIZ0N(0) = 150. cm.
MOISTURE CONTENTS
The terms soil moisture content and soil moisture movement are used
throughout this manual. More appropriately, the word "water" should
be substituted for "moisture." Although a matter of semantics, the
use of "water" is becoming more prevalent in literature.
Variables related to moisture contents include:
FCAP, card 1, group 1
BBC, TBC, TS, TD, SM, card 2, group 1
T0(I), group 4
similar variables in restart deck.
80
-------
All contents are on a volumetric basis, defined as:
9 - volume of water in sample
total volume of sample »• ••
Laboratory data are usually expressed gravimetrically as:
°f water in sample
rc;1
l J
dry weight of soil in sample
Both the volumetric and gravimetric moisture contents are
dimensionless quantities. They are related by the equation:
6 = aBD/Dw [6]
where
BD = bulk density of the soil (F/L3)*
Dw = density of water (F/L3)
In the cgs system, the density of water is approximately
1 gm/cm3. Expressing the bulk density in gm/cm3 yields:
8 = ctBD (numerically) [7]
The program uses volumetric moisture contents for input and output.
The above equations can be used to convert values. For example, if
the percent moisture by weight is 30 percent and the bulk density
is 1.5 gm/cm3, then:
9 = (0.30) (1.5) « 0.45 [gl
Since the bulk density of most soils is greater than 1, volumetric
contents will always be higher than the corresponding gravimetric
ones.
UNSATURATED FLOW PROPERTIES
The original program made use of statement functions and a series of
IF tests in the main program (5, p. 31) to compute the conductivity
and diffusivity. Since the IF tests appear in three places, making
changes in the program for different soils was not only a nuisance
* Dimensional Units are denoted in parentheses where L denotes length,
F denotes force, T denotes time, and - denotes a dimensionless quantity,
81
-------
but also an error prone process. To simplify program use and still
retain sufficient flexibility, the statement functions were removed
and the series of IF tests replaced by a single call to subroutine
PR0P. The general form of the call is:
CALL PR0P (ZK, TO, K, D)
where: 1 = volumetric moisture content at which the diffusivity
is computed
ZK = volumetric moisture content at which the conductivity
is computed
TD = volumetric moisture content below which the diffusivity
and conductivity are assumed negligible and set to
zero
K = conductivity at content ZK in cm/day (K is a real
variable
D = Diffusivity at content Z in cm2/day
The value of TD is obtained from the second card of group 1 and is
normally taken as the wilting point. Although the conductivity and
diffusivity are treated as constants during a time interval, they are
evaluated at different moisture contents. For flow between two nodes,
ZK is based on the average contents of the nodes at the start of the
time interval. Z is based on extrapolation of results from the last
two time levels for each of the nodes, averaged between them (5, equa-
tions [10], [11]). Under some conditions, it is possible for the
program to compute Z or ZK values greater than the saturation value,
TS on card 2 of group 1. Since TS is not passed to PR0P in the call
list, provisions should be included to limit the value of Z and ZK to
the saturated value when this occurs.
Subroutine PR0P can contain any degree of detail to approximate field
data. Methods may vary from using tabular values to employing exponen-
tial forms (5, p. 31). A sample is given in the section entitled
"Example of Subroutine PR0P."
The current form of M0ISTRE does not pass information to subroutine
PR0P specifying the depth or node number at which the properties are
computed. This is due to the fact that only homogeneous soils are
allowed. Original intentions were to allow for layered soils, hence
the structure of the group 3 cards. (As noted in 5, footnote A,
Table 23, p.33):
"Soil 'horizons' referred to in this context do not cor-
respond to morphological soil horizons. Program M0ISTRE
82
-------
models moisture movement in uniform soils only, so soil
'horizons' as referred to here have no function except to
correspond to the 'chemistry horizons' in the Biological-
Chemical Program for convenience."
The horizon references have been left in the program for future expan-
sion to handle variable flow properties. This would require the addi-
tion of the node number in the call list of PR0P. Properties could
then be related to a particular horizon.
A special form of PR0P has been developed for use with a specific
form of the moisture-release curve, as outlined by Campbell (23).
This method is purported to be as good as the Millington-Quirk
method (9) which is presently finding wide use and acceptance.
This routine receives some data through a C0MM0N block. The method
is detailed and is outlined as follows:
a. Plot moisture-release data on log-log paper.
b. If a straight line fits the plotted data reasonably well, then
the moisture-release curve is approximated by an exponential form
and the method of Campbell can be used.
c. Evaluate the slope of the line to obtain BEMP and the air-entry
potential AIRENT.
d. Select the option by setting IP0PT=1 on the first card of
group 1.
e. Include the group 8 card on input.
INITIAL BOUNDARY CONDITIONS
The top and bottom boundary conditions are established by variables
on card 2 of group 1. Values are set as follows:
a. Top boundary (TBC) -
if TBC = 1, set TBC = TM
if TBC * 2, set TBC = TS
if TBC = any other value, use TBC as read
b. Bottom boundary (BBC) -
if BBC - 1, set BBC = TM
if BBC » 2, set BBC = TS
if BBC = any other value, use BBC as read
83
-------
When initial conditions are not read from the group 4 cards or restart
data from cards or tape, all nodes other than the boundary nodes (i.e.,
internal nodes) are initialized to SM. Should either the group 4 cards
(CC=1) or restart data (IRESTR=1 or 2) be read, values must be included
for both boundaries. These override the TBC and BBC values at the
start of the first time step. When computations are finished for this
step, the top boundary condition depends on the availability of water
at the surface for infiltration while the bottom boundary condition
reverts to the BBC value (see MASS BALANCE). If water does remain at
the surface, the top boundary is set at TBC.
Logically, the group 4 input option (CC=1) and restart options (IRESTR=1
or 2) are not used for the same run. However, if they are, the initial
moisture contents read from the restart data override those from the
group 4 cards.
CONSUMPTIVE USE VALUES
Values read from the group 7 cards (CR0P=3) are the total consumptive
use values for each semimonthly period. The first period for a given
month extends from the 1st through the 15th while the second extends
from the 16th to the last day of the particular month. Consumptive
use is withdrawn at a uniform rate equal to the total divided by the
number of days in the given semimonthly period.
If the starting date of the run is not the 1st or the 16th, it is
necessary to adjust the input values to obtain the correct withdrawal.
Thus, if:
U = total required consumptive use for period (acre-inches/acre)
N = number of days in full period
A = actual number of days in period based on starting date
U = adjusted or input consumptive use (acre-inches/acre)
u = ^ m
For example, if the total consumptive use to be removed is 1.00 inch
and the starting date is July 13, then:
U = (LOO inch) (15) = 5>QO inches [1Q]
It is seen that the above equation merely converts the total consump-
tive use to the required withdrawal rate (U/D) and then computes the
input value by using the total number of days in the period. The
same approach is employed when the ending date of the run does not
fall on the 15th or the last day in the month.
84
-------
DEPLETED MOISTURE OPTION
A special option is available to handle winter precipitation for areas
experiencing an annual freeze-thaw cycle. It is selected by setting
IDEF=1 on the first card of group 1. The method is based on the
assumption that the entire profile is at least at a moisture content,
FCAP, normally considered to be the field capacity.
The starting day for a full year is considered as the thaw date while
the terminating day is the freeze date. No further infiltration or
redistribution of moisture occurs after the freeze date until the
thaw date of the next year. Winter precipitation during the period
is assumed sufficient to raise the moisture level throughout the pro-
file to the FCAP value. Any excess is presumed to melt and run off
or to evaporate.
The required quantity of water to raise the moisture content at each
node to FCAP is accumulated on the freeze date. If a node is already
higher, there is no reduction in the total. This total is printed on
output and used as the first application for the next year.
Water application data for the succeeding year must include the
depleted moisture as the first application. The total number of appli-
cations (APPS) must include it and the date and month specified for it.
The amount should be left blank. Although the application date
normally coincides with the thaw date, this is not necessary.
The depleted moisture option does not function until the end of the
first year's computations. Consequently, if the same technique is
employed for the starting year, the application must be computed and
included on the group 2 card. Similarly, restart data do not include
the quantity of depleted moisture from the previous year.
The following message is printed at the end of each year's computa-
tions when the depleted moisture option is used:
MOISTURE DEPLETED FOR YEAR = yy, MONTH = mm, DATE = dd (DAY NUMBER
= jj) IS = ii INCHES
where: yy = year number of run
mm = month number with mm=l for January, mm=2 for February,
etc.
dd = date within month
jj = Julian day number
ii = depleted moisture in the entire column in inches at
the end of the given date
yy, mm, and dd as well as jj correspond to the last day in a full year
or the stopping day of the run for the last year run.
85
-------
STARTING AND STOPPING DATES
Changes have been made to the program to allow runs beyond 1 year,
thus removing the limitations noted in Dutt, et al. (5, p. 35). All
dates are Julian with a perpetual calendar assumed (no leap years).
The starting year number of the run is given as YEAR (card type 2,
group 1) while the ending year is IST0P (card 1, group 1). These
years may be integer numbers starting with 1, 2, etc., or actual years
such as 1975, 1980, etc. Checks are made at the finish of each year's
computations to determine if the run is complete. Consequently, at
least a single year is run, even if IST0P is less than YEAR. In gen-
eral, the number of years in the run is (YEAR-IST0P+1).
The first year's computations commence with the Julian date LLSTRT
(card 1, group 1) and end with MM (card 2, group 1) unless YEAR =
IST0P. In this case the last day of computations is MMST0P (card
type 1, group 1). Subsequent years commence with date LL (card 2,
group 1) and end with day MM unless the last year (IST0P) is being
processed in which case the ending date is MMST0P.
As an example, if computations are desired to commence on July 12,
1975, and to end on November 14, 1977, with a thaw date of May 12 and
freeze date of November 27, the following values are used:
YEAR = 1975 IST0P = 1977
LLSTRT = 193 (i.e., July 12) MMST0P = 318 (i.e., Nov. 14)
LL = 132 (i.e., May 12) MM = 331 (i.e., Nov. 27)
The starting and ending dates for each year of the run are:
Year Starting date Ending date
1975 193 331
1976 132 331
1977 132 318
If there is no freeze-thaw cycle, then LL=1 (Jan. 1) and MM=365
(Dec. 31). Other variables are used in the manner outlined above.
DEFICIT MOISTURE
During a given time step, the program compares the quantity of water
specified for removal by the plants at a given node to the available
moisture. The available moisture is taken as the difference between
the moisture content at the start of the period and the input mini-
mum level, TD, usually taken as the permanent wilting point. Should
86
-------
the consumptive use exceed the available water, the difference is
computed and saved as the deficit moisture, and the consumptive use
actually removed is reduced to the available water. This technique
ignores the possibility of water moving from adjacent nodes during
the time step. Since contents are very low resulting in low trans-
mission rates, this appears to be a logical approach.
Should moisture levels rise above the minimum at later times, the
program again attempts to remove the specified consumptive use.
The method outlined above is obviously a simple approach to a complex
problem. The availability of water to a plant depends on soil mois-
ture levels, salinity conditions in the soil profile, crop type,
present health of the plants, stage of plant development, as well as
the physical soil properties. In addition, the relationship of the
entire root system to required water needs of the plant is involved.
When printed output is obtained at tenth-day intervals, a deficit
moisture condition results in the message:
DEFICIT MOISTURE SITUATION ENCOUNTERED nn TIMES THIS TENTH DAY.
AMOUNT IS aa CM
Where: nn is the number of individual situations that a
moisture deficit was encountered
aa is the accumulative deficit in cm3/cm2 or cm
A similar message appears with the daily output, except values corre-
spond to totals for the day and the message reads . . . TIMES THIS
DAY ... instead of ... TIMES THIS TENTH DAY ....
The cumulative deficit since the start of the run is DEFAMC which is
routinely printed on output, even if it is zero. It is also included
in the restart deck as is the cumulative deficit at a given node,
Z(J), since the start of the run.
MASS BALANCE
A mass balance check is included to assure program integrity and to
provide a means of assessing the adequacy of the time steps being
used in the simulation. The entire soil column is considered and an
accounting made of all inputs and outputs of water. Results of the
check appear with the printed output. The individual components are
carried forward in the restart data.
87
-------
The usual conservation of mass equation:
mass in - mass out = change in mass storage [11]
is applied to the soil column. The only water entering the column
is from infiltration of precipitation and irrigation water. Water
leaves the system as consumptive use and as deep percolation (leachate)
across the water table.
Defining:
ETS = cumulative infiltration (water percolating through the sur-
face into the column).
ET = cumulative plant evapotranspiration actually removed.
CL = cumulative water leached through the bottom of the column
(deep percolation). A positive value denotes water
leaving the column and a negative value denotes water
moving into the column from below the water table.
DIP = change in storage of water in the soil column since start
of run (final quantity less the initial quantity).
The units of the variables on output are in cm3 of water per cm2 of
surface area or simply cm.
The mass balance equation becomes:
ETS - CL - ET = DIP [12]
from which: CL = ETS - ET - DIP [13]
Letting: CHECK = ETS - ET - DIP [14]
CL = CHECK [15]
Printed output includes CL, CHECK, ETS, ET, and DIP. The CL and CHECK
values should compare reasonably well. Values for DIP and, consequently,
CHECK are only updated at the end of each day. Therefore when 0.1-day
output is selected, the mass balance check is only correct at the end
of each day.
The ET value is the quantity of water actually removed by the program.
An adjustment procedure is included to assure that this value equals
that on input when there are no deficits. This procedure is based on
the node spacing, root distribution with depth, and the program method
-------
for computing the withdrawal at each node. When low-moisture levels
occur and a deficit moisture (DEFAMC) results, the total consumptive
use on input should equal DEFAMC + ET. The mass balance equation is
unaffected.
A factor related to the mass balance computations is the approximate
manner of treating the upper boundary. When water remains at the
surface, the moisture content of the upper boundary is TBC. Once
all the water infiltrates, an extrapolation procedure is used to
estimate the top boundary moisture content for the next time step.
There is no assurance that the zero flow condition (i.e., a specified
gradient) is satisfied. Consequently, it is found that the program
computes water crossing the surface (in both directions) at later
times.
The ETS value is the total quantity of water crossing the surface and
is computed on the basis of gradients and soil properties as are all
flows. The quantity of water at the surface is stored in an accumu-
lator called HED which also appears in the restart data. When HED is
less than or equal to zero, all the water has infiltrated. However,
if additional water is computed as entering the column, it reduces
HED further (i.e., HED becomes more negative) and is included in both
ETS and a value called CI. Consequently, CI includes the net infil-
tration when water is at the surface plus any infiltration resulting
from errors in handling the upper boundary.
When there is no water at the surface and a quantity of water is com-
puted as leaving the column, flowing upward across the surface, it is
stored in a separate accumulator called CNA. It does not alter the
value of HED or CI, but will reduce ETS. The cumulative value of CNA
for the run is CNI.
From the above description:
ETS = CI - CNI [16]
It is emphasized that the CNA value represents the quantity of water
moving upward since all of the surface water from the last applica-
tion infiltrated the column. It is a result of the approximate treat-
ment of the upper boundary and as such is an error. The corresponding
error due to downward movement is the HED which increases negatively
until the next water application. When the application date is
reached, the program applies:
HED + CNA + AMT(I) I17!
HED and CNA tend to compensate. Together they represent the error
produced in handling the upper boundary which alters the application
-------
specified on input. After the last application is made for a given
year, the program checks the value of CNA. If it is appreciable,
(CNA greater than 0.05 cm), an application of HED + CNA is applied
the next day. This assures the value of ETS compares closely with
the accumulated AMT(I) specified on input. It is noted that units
on CNA, CNI, CI, and HED are all in cur of water per cm2 of area or
simply cm on the printed output.
RESTART OF PROBLEM
A problem may be run in successive stages using the restart capabil-
ity. Data necessary to restart a problem are automatically printed
during the run and may also be punched on cards or written to magnetic
tape (or disc storage) at the end of the run. The data can be used
in a subsequent run to continue the problem solution.
Restart data are directed to a permanent storage medium using the
IPUNCH option (card 1, group 1). If:
IPUNCH = 1 data are punched on cards
IPUNCH = 2 data are written to unit 10 which may be equipped to
a magnetic tape or to permanent disc storage
IPUNCH = any other value, no restart data are placed on a
permanent storage medium
Organizational structure on unit 10 and deck structure for the punch
file are shown in Figure 14. Tape 10 is unformatted (binary) while
the punch cards are formatted as indicated in Figure 14. The vari-
ables are defined below.
YEAR, M0NTH, IDTE The year, month, and date corresponding to
the time for which the restart data were
written. Month = 1 for January, 2 for
February, etc.
II The Julian day number corresponding to M0NTH
and IDTE.
IR Maximum rate of soil water movement between
nodes. It is used internally by the program
to select the time increment for solving the
finite difference equations (cm/day).
L Number of next irrigation or precipitation
event.
90
-------
B0T( beginning of tape) E0F(end of file mark)
| LOGICAL RECORD I | [LOGICAL RECORD 2 | | LOGICAL RECORD 3 | | LOGICAL RECORD 41 | LOGICAL RECORD 51
J L
J L
(TN(J),J = I,Q) (FN(J),J=I,Q) (ANT(J),J = I,Q) ( Z (J), J=I,Q)
YEAR, M0NTH,1DTE,II,CL,CHECK, IR,L,HED,C0NST,CI ,ETS,ET,CNA,CNI,DEFAMC
TAPE 10 RESTART DATA
(Z(J),J=I,Q)
6 VALUES PER CARD
4(Q-l)/6+ I Cords for each type of dateT/
/_
(ANT(J),J = I,C
/
/(FN(J),J=I,Q)
rN(J),J = l,Q)
/
DEFAMC
C0NST,C!,ETS, ET, CN A.CNI
/YEAR,M0NTH,IDTE ,CL,
CHECK,I R,L,HEAD
FORMAT (6EI3.6)
FORMAT (15,213,14, 3E,13.6,13,EI3.6)
FIGURE 14. PUNCHED CARD RESTART DECK-UNSATURATED FLOW SUBMODEL
( RESTART DATA)
91
-------
C0NST
Q
TN(J)
FN(J)
ANT(J)
Z(J)
Initial quantity of water in the unsaturated
soil column (cm3 of water per cm2 of surface
area or cm).
Number of nodes in the unsaturated column used
to solve the finite difference equations.
Present moisture content at node J
(dimensionless).
Present rate of water movement between nodes
J and J + 1 (i.e., between node J and the next
one below it) (cm/day).
Anticipated moisture content for the next time
step to be used in evaluating the diffusivity
at node J (dimensionless) (see UNSATURATED
FLOW PROPERTIES).
Cumulative deficit moisture at node J (cm3
of water per cm2 of surface area or cm).
CL, CHECK, HED, CI, ETS, ET, CNA, CNI - See MASS BALANCE.
DEFAMC See DEFICIT MOISTURE
Restart data are automatically written to the common output device
(printed output) regardless of the value of IPUNCH. Data are written
at the end of 15th day, at the end of each month, and at the end of
the run, regardless of the termination date. This output may be used
to generate the required restart deck if an abnormal termination
occurs. Items corresponding to those in logical records 2-5 are
listed in nodal order, starting at the surface with seven values per
line.
Restart data generated by the IPUNCH option are used by the program
by setting the restart flag IRESTR appropriately. If:
IRESTR = 1 data are read from cards (Figures 13 and 14)
IRESTR = 2 data are read from unit 10
IRESTR = any other value, restart data are not used.
Restart data are intended to be used for the day following the date
on which it was written. The program does not check dates although
they do appear with the printed input data.
When a problem is restarted and tape 5 output is saved, it is necessary
to position tape 5 correctly (see OUTPUT).
92
-------
When restart data are read, the CC=1 option should not be used. How-
ever, if it is, moisture contents from the group 4 cards will be over-
ridden by the restart data (see INITIAL AND BOUNDARY CONDITIONS). Run
information and options (group 1), water applications (group 2), hori-
zon data (group 3) as well as root zone data (group 6), consumptive
use values (group 7), and unsaturated flow properties (group 8) when
used must be included. Starting and ending dates as well as water
applications for the entire year must be included, corresponding to
the year in which the restart data were developed.
LIMITATIONS
The program has been developed with the following restrictions:
a. Only homogeneous soils are allowed. The extension to layered
conditions would be relatively simple (see UNSATURATED FLOW
PROPERTIES).
b. The bottom boundary is at a fixed location. The fluctuating
water table is not treated.
c. A single, constant root pattern is used.
Problem size is restricted by allocated storage to:
a. maximum number of nodes = 60
b. maximum number of water applications = 50
c. maximum number of horizons = 9
There is no limit on the number of years run.
Use of the CR0P = 1 or 2 option is not recommended since the Blaney-
Criddle coefficients are for two crops based on Arizona data (see
DATA DECK STRUCTURE).
OUTPUT
Program output consists of printed listings, punched cards, and both
computed results and restart data written to magnetic tape (or disc
storage).
a. Printed data are listed by the program on the file called
0UTPUT. It includes both standard and optional output.
93
-------
(1) Standard output includes:
(a) Input data: unsaturated flow properties when Campbell's
method is used, water application dates and amounts, evapo-
transpiration values, and restart data read from tape or
cards when the run is restarted.
(b) Computed results: the amount of water applied on a given
date; mass balance results, upper boundary information, def-
icit moisture values and the number of time increments used
for a given day for each day of the run; and restart data at
the end of the 15th, the month, and the run.
(2) Optional printed output includes:
(a) Input data: control and run parameters, initial mois-
ture contents, soil horizon data, and root distributions if
AA = 1 on card 2, group 1.
(b) Computed results: moisture contents on a daily basis if
BB = 1 on card 2, group 1; mass balance results, upper bound-
ary information, deficit moisture values, the number of time
increments used thus far for the day, and moisture contents at
0.1-day intervals if ITENTH = 1 on card 1, group 1.
(3) Special printed messages include:
(a) Depleted moisture message at the end of each year and at
the end of the run when IDEF = 1 on card 1, group 1 (see
DEPLETED MOISTURE OPTION).
(b) Deficit moisture message whenever a deficit is encoun-
tered during a given day (see DEFICIT MOISTURE).
(c) Unstable solution message whenever this situation is
encountered (see HINTS ON USE).
(d) Error messages regarding tape 5 when a saved tape is
used (ISAVE = 1 on the group 1 card) and the error condition
is detected (see discussion below on tape 5 output).
The deficit moisture and unstable solution messages are given
at the end of each 0.1 day when ITENTH = 1 and the given situ-
ation occurs.
b. Restart data are included in the printed output and may be
written to unit 10 or to the punch file (see RESTART OF PROBLEM).
94
-------
c. Computed results are written by the program on unit 5 as stand-
ard output at the end of each 0.1 day. Data are written without
format in binary form by the single write statement:
WRITE (5) YEAR, II, XT, CI, CL, HED, ETS, DEPAMC, (J, TN(J),
Z(J), SF(J), U(J), J=l, Q)
Where: YEAR = current year number
II = Julian day number
XT = time within the day corresponding to the
written data (i.e., 0.1 day, 0.2 day, . . .
1.0 day)
CI, CL, HED, ETS - See MASS BALANCE
J = node number (see NODE SPACING)
TN(J) = volumetric moisture content at time XT for
node J (cm3/cm3 or -)
Z(J) = deficit moisture at node J during the last
0.1-day interval (cm3 of water/cm2 of sur-
face area or cm)
SF(J) = total soil moisture flux (i.e., total quantity
of water moving) between node J and (J+l)
during the last 0.1 day (cm3 of water/cm2 of
surface area or cm)
U(J) = evapotranspiration rate for node J based on the
input data (cm3 of water/cm2 of surface
area/day or cm/day)
An end-of-file mark is written at the end of each year's data. Con-
sequently, each year's data can be considered a single logical file.
When a problem is restarted and a saved tape 5 is to be written to,
it is necessary to position it for the initial write.
Use of a saved tape is indicated by setting ISAVE=1 on the group 1
card. The initial file to which additional computed results are to
be added is denoted as INFIL5 on the group 1 card. The program skips
over the first (INFIL5-1) files, positioning tape 5 to the first
record in INFIL5. If the starting date LLSTRT equals LL (the first
day in a full year), the tape is assumed positioned and ready to go
(i.e., INFIL5 is empty). If LLSTRT is not equal to LL, 10 records
per day are read, the date read is checked against (LLSTRT-1), with
the process stopping when it is found. The next write will then
generate data for the first 0.1-day interval in day LLSTRT.
Two abnormal situations can occur during the process, resulting in
error messages and program termination:
a. The end-of-file is reached before day number (LLSTRT-1) is
reached. Restarting at LLSTRT would leave a gap in the records
95
-------
in INFIL5. This situation includes the case of an empty file.
The following message is printed:
END FILE FOUND ON TAPE 5 BEFORE DAY NO. LLSTRT-1 FOUND,
EXECUTION TERMINATED
b. The first day number read from the file may exceed LLSTRT.
It is then impossible to position tape 5 correctly. The follow-
ing message is printed:
DAY READ FROM TAPE 5 EQUALS OR IS GREATER THAN STARTING
DAY, EXECUTION TERMINATED
Tape 5 is intended to be used as input to an interface program which
interfaces the unsaturated flow program with an unsaturated chemistry
and a drainout program.
HINTS ON USE
The following items are intended to guide the user in the intelli-
gent use of the program. Some of the items are discussed elsewhere
in the user's manual but are noted here for ready reference.
Selection of Maximum Time Step
Mass balance errors generally increase as the average time step for
solving the finite difference equations increases. This is due to
errors inherent in using finite differences as well as those result-
ing from the techniques that are used to handle the nonlinear soil
properties.
Minimum computer costs dictate as large a time step as practical.
The time step is automatically reduced when water is applied. This
reduction is required because moisture conditions and, consequently,
the unsaturated flow properties are rapidly changing with time. The
reduced step is 0.001 day corresponding to 1,000 steps per day. These
steps are gradually increased as the maximum rate of moisture movement
in the profile decreases. The user limits the maximum size to M steps
per day on card 2, group 1.
Since results are written to tape 5 on a 0.1-day basis, the largest
time step is logically 0.1 day (M=10). However, the program checks
the size of the step it is about to use and reduces it, if necessary,
to obtain computed results at the 0.1-day interval.
Because of this, situations can occur in which successive time steps
exhibit large variations in size. For example, steps of 0.045, 0.050,
and 0.005 might result from specifying M=20. Experience has shown that
96
-------
better results are obtained when the tenth-day period is divided into
more nearly equal steps, e.g., 0.035, 0.035, and 0.030. This is due
to the method of handling the upper boundary and of estimating the
anticipated moisture content used in evaluating the diffusivity f5
equations [10], [11], p. 7).
For most soils, a value of M=30 is recommended. The actual number
of time steps used is given as NSTEPS in the printed output. Both
mass balance results and conditions at the upper boundary should be
used in selecting the appropriate value of M.
Unstable Solution
Because of the methods used to deal with the nonlinear flow proper-
ties, the approximate treatment of the upper boundary, and numerical
solution errors, the computed moisture contents may exceed the physi-
cal bounds. When the computed value is less than zero or greater than
saturation (i.e., TS) at a node for one time step, a counter is incre-
mented. The total number of such occurrences is noted in the printed
output for each 0.1 day when ITENTH=1 and for each day. The printed
message is:
UNSTABLE SOLUTION SITUATION ENCOUNTERED nn TIMES THIS TENTH
DAY or (. . . THIS DAY)
Normally, this situation only occurs on application dates when the
soil profile is very dry and the amount of applied water large. Mois-
ture contents then exceed TS. Under rare conditions, moisture levels
may fall below zero, especially when the wilting point, TD, is very
near zero. In either case, if the count is low and computed contents
exceed the limits only slightly (say 0.01 or less), the situation is
not serious and can be ignored. Generally, major difficulties result
only when very heavy soils (low conductivity) are encountered. In
this case, improvement can be made by increasing the number of nodes
used to simulate the soil column.
Adequacy of Unsaturated Flow Properties
The unsaturated flow properties should be evaluated by performing sev-
eral simulations to duplicate soil properties commonly used by soil
scientists. These include field capacity, specific yield, wilting
point, and water holding capacity.
Uniform Initial Conditions
Whenever uniform initial conditions are desired, the value of CC on
the type 2, group 1 card should be left blank or set to zero. The
97
-------
SM value on the same card is then used for initial conditions at all
internal nodes.
Deficit Moisture
Whenever the deficit moisture situation is encountered, the reason for
it should be determined. Normal irrigation practices attempt to pre-
vent such conditions. Items that should be considered include the
unsaturated flow properties, consumptive use rates, and the irrigation
schedules.
Boundary Conditions
The top and bottom boundary conditions need not be set at saturation.
The top boundary condition, used when water is at the surface, can be
reduced to account for surface resistance due to compaction, inwashing,
swelling of colloids, and the breakdown of soil structure. The bottom
boundary may not represent a water table. Consequently, any moisture
condition can be used to model it, corresponding to the desired tension.
EXAMPLE OF SUBROUTINE PR0P
Introduction
When the method of Campbell is not used to compute the unsaturated flow
properties, the user must supply his own property subroutine PR0P. The
general form of the call is:
CALL PR0P (Z, ZK, TO, K, D)
where: Z = volumetric moisture content at which the diffusivity
is computed
ZK = volumetric moisture content at which the conductivity
is computed
TO = volumetric moisture content below which the diffusivity
and conductivity are assumed negligible and set to
zero
K = hydraulic conductivity at content ZK in cm/day (a real
variable)
D = diffusivity at content Z in cm2/day
Subroutine PR0P may use any method to determine the conductivity and
diffusivity. Tabular values with an interpolation function or mathe-
matical relationships may be used. Particularly convenient mathemati-
cal relations are those based on exponential forms.
98
-------
Exponential Forms: Semilog
Consider the case where data plotted on semilog paper can be repre-
sented by a straight line as shown in Figure 15. Points 1 and 2 are
arbitrary points on the straight line while r is a reference point at
which the unsaturated flow property is Pr. The equation of the
straight line is:
log P = log Pr + slope X (6-6r) [18]
The slope is given by:
log P2 - log P! log (P2/PX)
slope = - - - - - — — [19]
(62-61)
Equation [18] becomes:
P . log pr . (e-er) [20]
The conversion from base 10 to base e logarithms is:
Iog10 x = 0.43429 In x or Inx = 2.303 Iog10x
Consequently, equation [20] becomes:
ln /P \. 2.303 l!!JV!ll (e-8r) [21]
The logarithmic function x=logby is defined as the inverse of the
exponential function and is given by the relationship y=b* where
b is the base. Thus for natural logarithms b=e, x=lny, and y=ex.
Consequently, equation [21] becomes:
[2.303 log (P2/Pi) (e-6r)]
- -—
The exponent of e can be written as:
E « C (8-er)
where: C = 2.303 2—±- t24!
99
-------
>
^-
-------
C is a constant (essentially the slope) which is evaluated by select-
ing two points 1 and 2 so that 62>81 and obtaining the 6's and corre-
sponding P values. In addition, a reference point is selected (which
could coincide with either point 1 or 2 if desired) and the values of
6r and Pr obtained.
Exponential Forms: Log- log
Consider the case where data plotted on log-log paper can be repre-
sented by a straight line. The equation of the straight line is:
fp i !og (P,/?,)
lnl ' 2-303 log (e/v [25'
[2.303 log (P2/Pi) ]
P=Pr e - - 2 IJ log (0/er) [26]
log (62/ej) r
The exponent of e can be written as:
E = c log (e/er) [27]
where: C = 2.303 1<>g ^f1^ [28]
log (Oj/ep J
C is evaluated in a similar fashion to that for the semilog form.
Example
The log-log form is often preferable to the semilog form as it usually
results in a better straight line fit. Several straight line segments
may be required to approximate actual data, regardless of the form.
The following example illustrates use of the semilog form with several
straight line segments.
Measured data defining the moisture-release curve (h vs. 6) and the
conductivity (K vs. 0) are plotted as points on the semilog plot of
Figure 16. Possible functional relations are shown by the solid
curves. Diffusivity values were computed by selecting a given 6,
determining K(6), evaluating the slope of the moisture-release curve,
and applying the relationship:
D(6) = K(6) [29]
101
-------
,000,000
500,000
.002
.001
0 O.I 0.2 0.3 0.3"
0 VOLUMETRIC MOISTURE CONTENT (-)
FIGURE 16. UNSATURATED FLOW PROPERTIES
102
-------
Portions of the moisture-release curve exhibit very steep slopes making
it difficult to evaluate 9h/36. At best, the procedure is approximate?
The dashed curves represent poor, though possible, approximations of
the functions. They are used to illustrate the procedure. A single
reference point is chosen for the conductivity functions:
6r = 0.166 where K,. = 1.2 cm/day
Different reference points are chosen for each segment of the diffu-
sivity relationship:
6r = 0.166 where Dr = 158.75 cm2/day
6r = 0.255 where Dr = 550 cm2/day
In general, reference points are chosen at the intersection of two
straight line segments. This assures that values picked from the
curve will produce the same value at the intersection by either rela-
tionship when the equations are used. Although this is good practice,
it is not absolutely necessary.
In the Fortran coding for subrouting PR0P, it should be noted that the
property relations can be extended below the TD value. In this case,
the TD value passed in the call list is simply not used.
CAMPBELL'S METHOD FOR DETERMINING UNSATURATED FLOW PROPERTIES
Introduction
Gaylon S. Campbell(23) has proposed a simple method for developing the
functional relationship between the moisture content and hydraulic con-
ductivity. It is based on a relatively simple form of the moisture-
release curve. Campbell found that agreement between conductivities
calculated by this method with those experimentally determined for five
soils was as good as with the Millington-Quirk method.
The soils tested included soils classified as loams, sandy loams, and
sands. Saturated conductivities ranged from 0.0012 cm/min to 1.12 cm/
min. Saturated volumetric moisture contents ranged from 0.35 to 0.52.
Because this method is simpler to use than the better known and widely
accepted Millington-Quirk method, its use is appealing. In addition,
it is easily extended to include the diffusivity function.
Theory
If a soil sample is saturated with water, the corresponding soil-water
potential is atmospheric. As water is removed from the sample, the
103
-------
potential falls below atmospheric. Using atmospheric pressure as
the gage reference, the soil water potential is negative. This is
commonly denoted by referring to them as tensions.
A plot of the soil-water potential as a function of the volumetric
moisture content is variously known as the moisture-release curve,
moisture retention curve, or the moisture characteristic. If it
can be represented by the functional relationship:
3r
\
the hydraulic conductivity is:
2b+3
KQ = K- [31]
where: ty is the soil-water potential (L)
$Q denotes ij> at the arbitrary content 0 (L)
i|»r denotes \J» at the reference level 6r (L)
6 is the volumetric moisture content (L3/L3 or -)
6r denotes 0 at the reference level r (L3/L3 or -)
KQ is the hydraulic conductivity at
content 0 (L/T)
Kj, is the hydraulic conductivity at
the reference content 0_ (L/T)
b is an empirical constant determined
from the data (-)
Equation [31] is based on the following assumptions:
a. Flow of water in the soil is controlled by the smallest of two
pores in sequence.
b. Only pores that are connected together contribute to the total
hydraulic conductivity (i.e., isolated and dead-end pores are not
effective.
c. Pores fit together randomly.
Based on these assumptions, the pore radii and pore size distribu-
tion function are used to determine the conductivity. The pore radii
are related to water contents by the moisture-release curve and the
capillary-rise equation. A pore interaction term is used to modify
the final conductivity function.
104
-------
In addition to the conductivity, solution of the unsaturated flow
equation requires the diffusivity. It is defined as the conductivity
at a given moisture content divided by the corresponding slope of the
moisture-release curve:
[32]
where: DQ is the diffusivity at content 6 (L2/T)
Based on equations [30] - [32], the diffusivity becomes:
[33]
6r
The diffusivity at the reference level r can be found from equa-
tion [33] by setting 6 = 6r. Then:
Dr " -T-^ [34]
r
Dimensionless Forms
The preceding equations are easily recast into dimensionless form
as follows:
£• (tr
c-fS-)**'
^ IV
n /« \ k*2
£L. /1_\ [37]
Dr ^ry
They are useful in comparing properties of different samples. Plots
of each dimensionless property on log-log paper may be inspected for
similarities and differences.
105
-------
Evaluation of Empirical Constant b
It can be shown that a relationship of the form given by equation [30]
will plot as a straight line on log-log paper. The slope is equal to
the value of the empirical constant b. The following procedure can be
used to determine b from a set of moisture retention data.
a. Plot ^ versus 0 on log-log paper.
b. Fit a straight line through the data points.
c. Evaluate the fit. If it is unacceptable, Campbell's method can-
not be used. An alternate method such as Millington-Quirk should be
tried.
d. Select two points on the straight line. (Point 1 corresponds to
the highest 6.)
e. Calculate b by either of the following equations:
log OP o/if*,)
or
b = -
log (e2/ei)
In
In
[38]
[39]
Note that the base of the logarithm is irrelevant and b should
be a positive quantity.
The 6's in the previous equations are volumetric. In general, labora-
tory data are reported on a gravimetric basis.
Thus:
volumetric - 6 = volume of soil water/total volume of sample
(L3/L3 or -)
gravimetric - a = weight of soil water/dry weight of soil
(F/F or -)
It can be shown that:
9 = a BD/D
w
[40]
106
-------
where: BD = bulk density of soil (i.e., dry weight of soil/unit
volume of total sample) (F/L3)
Dw = density of water (i.e., weight of water/unit volume
of water) (F/L3)
Substitution of equation [40] into equation [38] yields:
log
b =
log
Consequently, determination of b can be based on either volumetric or
gravimetric data. Similarly, any set of consistent units for the poten-
tials may be used such as centimeters, inches, feet, or bars* since con-
version from one set of units is merely ip multiplied by a constant in
which constants in numerator and denominator cancel.
Significance of Reference Level r
A specific state, termed the reference level, is noted and used in many
of the previous equations. Certain variables are referenced to this
level, as denoted by the subscript r. The reference level is merely
the moisture level at which the hydraulic conductivity is known. Prac-
tically, the saturated conductivity K is most easily measured in the
field or laboratory. Then
Kr = Ks [42]
6r = 9S [43]
Since Dr is computed from the above variables, it can be excluded. 9S
is the moisture content at saturation and tye is the corresponding poten-
tial. ipe is called the air-entry potential. Its physical significance
and determination are discussed in the next section.
Air-entry Potential, $Q
If a soil is saturated and at equilibrium with free water at the same
elevation, the pressure is atmospheric and the tension head is zero.
* 1 bar = 1 atmosphere 14.7 lb/in2 which is equivalent to a column of
water = 33.92 feet = 407.1 inches = 1,033 cm.
107
-------
Hillel (24) states:
"If a slight suction, i.e. a water pressure slightly
subatmospheric, is applied to water in a saturated soil,
no outflow may occur until, as suction is increased, a
certain critical value is exceeded at which the largest
pore of entry begins to empty. This critical suction
is called the air-entry suction. Its value is small in
coarse -textured and in well -aggregated soils. However,
since in coarse- textured soils the pores are often more
nearly uniform in size, these soils may exhibit critical
air-entry phenomena more distinctly and sharply than do
fine-textured soils."
Carefully run laboratory tests on undisturbed samples might provide
a value for i{)e. This has never been done by the Bureau of Reclama-
tion. However, a suitable value might be obtained by extending the
straight line of the log-log plot of potential versus moisture con-
tent to saturation.
A review of current literature reveals a lack of numeric values for
the air-entry potential. However, values should range from 5 to
100 cm depending on the soil texture and pore size distribution.
Lacking field data, the following values are suggested:
Texture Range in i|>e
sand 1- 5 cm
silt 10-20 cm
clay 40-60 cm
It is noted that if>e cannot assume a zero value since the relationship
of equation [30] would be meaningless.
The air-entry potential concept is related to capillary forces. This
suggests the use of the capillary rise equation for cylindrical tubes
for soils on which a mechanical analysis is run. Based on a water tem-
perature of 25° C (77° F) , an assumed contact angle of zero, and the
previously suggested air-entry values, the following equations are
obtained :
reff = (10PS +3.3 PSI + Pc)/34,000 [45]
ipe - 0.147/reff [46]
108
-------
where: reff = effective capillary radius of pores (cm)
PS> PSI> PC = percent of sand, silt, and clay, respectively,
based on the mechanical analysis
ij)e = air-entry potential (cm)
An effective capillary tube radius is first computed using equa-
tion [45]. Since the larger pores drain first, this equation places
the most weight on the sands and the least weight on the clays. The
capillary rise equation [46] is then used to compute the air-entry
value. If the soils are of uniform texture, the following values are
given by the equations:
Texture reff ^e
(100%) (cm) (cm)
sand 0.0294 5
silt 0.00971 15
clay 0.00294 50
There is no provision in equations [45] and [46] to include soil struc-
ture. At best, their use provides approximate air-entry values when
measured data are lacking.
109
-------
SECTION XI
USER'S MANUAL - DRAINOUT PROGRAM
INTRODUCTION
This program was designed to predict the response of a subsurface
drainage system consisting of parallel, equally spaced tile drains
to deep percolation inputs. Both the water table elevation midway
between drains and the drain discharge are computed. In particular,
this program was developed as part of a water quality simulation
model with the computed monthly drain discharge acting as input to
other programs.
The program may be used in two distinct ways:
a. Analysis. - An existing or designed drainage system can be
analyzed for its response to a given set of deep percolation inputs.
b. Design. - Given the porous media properties, depth to barrier,
effective drain radius, depth to drains, and the deep percolation
schedule, the drain spacing that will adequately provide the speci-
fied drainage depth can be determined.
Techniques conform closely to those used routinely by the USER in deter-
mining drain spacings. Details are given in the program description.
Deep percolation data can be inputed by either of two methods:
a. Magnetic tape input of daily values based on the results of
an unsaturated flow model. This tape is written by an interface
program.
b. Card input of values based on any desired method. They may
be developed by considering crop type, agricultural practices,
climatic conditions, and field experience.
Output includes listed results, an optional magnetic tape for use
with a cathode ray tube plotting routine and optional punched out-
put for use with the drain effluent prediction portion of the water
quality model. (CRT plotting routines are not included in this
report.)
UNIT NUMBERS
Five input-output devices are referenced by the program:
110
-------
Unit Unit description
designation or purpose
1 The unit on which the magnetic tape containing deep per-
colation data is mounted. These data are written by the
interface program and computed by the Unsaturated Flow
Program.
3 The unit to which computed results are written on a
saved magnetic tape. An entirely new tape may be writ-
ten or results for additional years added to a previously
generated tape. It is used in conjunction with the util-
ity program to obtain plots of the results.
5 The common input tape or unit (CIT). It is the computer
system's standard input device.
6 The common output tape or unit (C0T) 4 ];•£ ^s tfte com-
puter system's standard output device.
7 The punched card output device of the computer system.
DATA DECK STRUCTURE
There are six groups of input cards in the data deck, as shown in
Figure 17.
Cards in groups 1 through 5 apply to both analysis and design. Most of
the variables on these cards are identical for either use. However,
certain variables are used in slightly different ways, others have
unique restrictions, and some are not used for design. These are pref-
aced by an asterisk. Their use in design is discussed separately along
with the group 6 card in the section, Design Inputs.
CARD GROUP 1 - Title Card
Column Var Var Variable
No. type name description
Card 1 of 1
1_80 A TITLE An 80-character alphanumeric title used
to identify all printed and plotted out-
put. If this title is desired to be
centered on output, it must be centered
on the card.
Ill
-------
Repeat group 5 for each year
of deep percolation data.
Years deep percolation
/I Years deep percolation
Run options and
length
Porous media properties
Drain geometry and
water table data
Title card
FIGURE 17. DATA DECK STRUCTURE - DRAINOUT SUBMODEL
112
-------
Column
No.
Variable
description
CARD GROUP 2 - Drain Geometry and Water Table Data
For further clarification, see Figure 1. The water table position
referred to by YSTRT, YEST, YLIM, and YTEST means the height of the
water table above the horizontal centerline of the drains at a point
midway between adjacent drains.
Card 1 of 2
1-10
11-20
21-30
31-40
41-50
51-60
61-70
71-80
R SPACE Drain spacing (feet).
R D Distance from horizontal centerline of
drain to barrier (feet).
R ER Effective radius of drain (feet).
R CLSTRT Amount of deep percolation added to the
aquifer on the first day (January 1)
of the initial year of the run (acre-
inches/acre or inches). (See RESTART
OF PROBLEM.)
R YSTRT The water table position at the start of
the initial year of the run (feet).
This value can be negative, zero, or
positive. (See NEGATIVE STARTING POSI-
TION, RESTART OF PROBLEM.)
R YEST The estimated maximum water table posi-
tion (feet). (See SELECTION OF DRAIN-
AGE CASE.)
R YLIM The water table position below which the
drainout is considered complete (feet).
When the computed position is below this
value, it is set to zero for output as
is the discharge until another deep per-
colation occurs or the run ends.
R YTEST A test value used to determine when
dynamic equilibrium is reached (feet).
If the water table is within YTEST at
the end of successive years, dynamic
equilibrium is assumed. (See DYNAMIC
EQUILIBRIUM.)
113
-------
Column
No.
Card 2 of 2
1-10
11-20
21-25
Var Var Variable
type name description
R DCSTRT The value of the cumulative leachate or
deep percolation at the start of the
initial year when input is via magnetic
tape 1 (cm3/cm2 or cm). (See DEEP
PERCOLATION DATA.)
R YBLIM The absolute value of an instantaneous
buildup (deep percolation/specific
yield) which must be equaled or
exceeded before it is applied (feet).
(See APPLICATION OF METHOD.)
I MXTDR The maximum length allowed for a drainout
period (days). (See APPLICATION OF
METHOD.)
CARD GROUP 5 - Porous Media Properties
Card 1 of 1
1-10
11-20
R PERM Aquifer permeability (feet/day).
R SY The storage coefficient expressed as a
decimal (ft3/ft3 or dimensionless). It
is often taken as the specific yield.
It is the ratio of drainable or fill-
able voids volume to the total volume.
CARD GROUP 4 _-_ Run Opt ions and Length
Card 1 of 1
1- 5
NYRS The number of years of
data to be processed
Each year is one file
input is by magnetic
of group 5 cards when
(See ICRDIN below.)
NYRS may be zero when
is used to calculate
table.
deep percolation
for this run.
on unit 1 when
tape or one set
input is by cards.
It is noted that
the N0PT option
the falling water
114
-------
Column
No.
6-10
Var Var
type name
I INFILE
11-15
I N0PT
16-20
I NREPT
21-25
I NSAVE
26-30
I NTYPE
Variable
description
The number of the initial file to be
read on unit 1 when deep percolation
inputs are by magnetic tape. (See
RESTART OF PROBLEM.)
Option to allow calculating drainout
after the last year of deep percola-
tion input data.
If N0PT ^ 1, the program stops at the
end of the last year
having inputs.
If N0PT = 1, the program continues to
compute results using
zero inputs until drain-
out is completed accord-
ing to YLIM.
Option to allow repetition of the last
year's deep percolation input data for
successive years. The option is avail-
able for both card and tape input.
Note that N0PT may be used in conjunc-
tion with NREPT.
If NREPT *-
not used.
If NREPT > 0, drainout calculations
are continued with the last year of
input repeated NREPT times, unless
dynamic equilibrium is reached. (See
DYNAMIC EQUILIBRIUM.)
Option to allow writing of daily and
monthly results on unit 3.
If NSAVE = 1, data are written. (See
OUTPUT.)
Selects the drainage case to be used in
the analysis.
If NTYPE = 1, analyze as drains above
barrier.
115
-------
Column Var Var Variable
No. type name description
If NTYPE = 2, analyze as drains on
barrier.
If NYTYPE = any other value, use the
program criteria to
select the case.
(See SELECTION OF
DRAINAGE CASE.)
31-35 I IYEAR The year number corresponding to the
initial year being processed in the
current run. It may be 1, 7, 1972,
etc., and is used to identify output.
36-40 I NPUN Option to allow punching of the average
monthly discharge for each year of the
run.
If NPUN = 1, then punch cards. (See
OUTPUT.)
41-45 I INFIL3 The number of the initial file to be
written on unit 3. (See RESTART OF
PROBLEM.)
46-50 I ICRDIN Option to allow input of deep percolation
data by cards.
If ICRDIN = 1, then input is by cards,
group 5.
51-55 I IDESGN Option to specify the design mode.
If IDESGN = 1, the program is to be used
for design.
If IDESGN ± 1, the program is to be used
for analysis.
CARD GROUP 5 - Deep Percolation Data (ICRDIN = 1 Only)
Card 1 of ?
11-12
IM0N Number corresponding to month of appli-
cation (i.e., January is 1, June is 6,
etc.).
116
-------
Column Var Var Variable
No- type name description
13-14 * IDAT The date of the month of application
(i.e., an application on February 12
is IM0N = 2, IDAT = 12).
15 A ADENT A single alphanumeric character which can
be used to designate the type or source
of deep percolation (i.e., I for irriga-
tion, P or R for precipitation, etc.)-
It is printed on output.
16-20 R AMT Amount of deep percolation on the given
date (acre-inches/acre or inches).
Note: The IM0N, IDAT, ADENT, AMT is repeated for each application
with seven values per card. The applications need not be in chrono-
logical order, except they must be in groups by years, with the first
application being first and the last application being last for each
year. The end of a year's input is indicated by a blank IM0N, IDAT,
ADENT, AMT set. If the number of applications in a year is a multiple
of 7, then an additional blank group 5 card is required. (See DEEP
PERCOLATION DATA.) A perpetual calendar is used so that no application
should be made on February 29.
The deep percolation inputs for a given day are assumed to occur
throughout the day. They are applied as an instantaneous deep perco-
lation at 12 midnight (i.e., at the start of the next day). Conse-
quently, an initial application on the first day of the initial year
must be specified by the value of CLSTRT (card group 2).
DESIGN INPUTS
The section entitled Design Method should be consulted for clarifica-
tion of the design inputs.
CARD GROUP 2
Card 1 of 2
1_10 R SPACE Initial trial drain spacing (feet). If
the Donnan formula is used, SPACE should
be left blank.
117
-------
Column
No.
51-60
CARD GROUP 4
Card 1 of 1
1-5
R
VEST
NYRS
11-15
16-20
N0PT
NREPT
21-25
31-35
36-40
NSAVE
IYEAR
NPUN
Variable
description
The design maximum water table position
(feet).
The number of years of deep percolation
data to be processed for this run. For
design, only a single year of inputs is
allowed. Therefore, NYRS must equal 1.
(See INPUT CHECKS AND ERROR MESSAGES.)
This option has no significance when the
design mode is being used and is ignored.
The maximum number of years allowed to
reach dynamic equilibrium for one trial
spacing. (See INPUT CHECKS AND ERROR
MESSAGES.)
For design, only the results for the cor-
rect drain spacing (i.e., the designed
spacing) at dynamic equilibrium are
available for output. If the spacing
is not determined, there are no water
table or drain discharge results.
Therefore:
NSAVE. - The daily and monthly
results at dynamic equilibrium
for the designed system are writ-
ten on unit 3 when NSAVE = 1.
IYEAR. - The year number associated
with the results at dynamic
equilibrium.
NPUN. - The average monthly discharge
values for the designed system at
dynamic equilibrium are punched when
NPUN = 1.
118
-------
Column
No.
41-45
Var
type
Var
name
INFIL3
Variable
description
The number of the initial file to be
written on unit 3. It is noted that
multiple files may be written on unit 3,
even though each file represents a dif-
ferent drainage system.
CARD GROUP 6 - Design Data (IDESGN = 1 only)
Card 1 of 1
1- 5
6-10
11-15
16-20
21-25
26-35
I
R
36-45
R
IM0NST Number of month corresponding to start
of period from which SEEP is computed
(i.e., January is 1, June is 6, etc.).
IDATST The date corresponding to the start of
period given by IM0NST.
IM0NFI Number of month corresponding to end of
period from which SEEP is computed.
IDATFI The date corresponding to the end of
period given by IM0NFI.
NTRYMX Maximum number of trial spacings allowed.
YCHECK The absolute value of the difference
between the computed maximum rise at
dynamic equilibrium and the design
value YEST which cannot be exceeded
when the drains are properly spaced.
SEEP Steady-state seepage rate to be used
when the Donnan formula is used for the
initial spacing (feet/day).
SELECTION OF DRAINAGE CASE
The USBR drainage design procedures make use of the two limiting cases:
a. Drains an infinitely long way above an impervious barrier.
b. Drains resting on an impervious barrier.
119
-------
Any real system must fall between these extremes. Based on experi-
ence, the USER has found the following criteria useful in selecting
the appropriate case:
a. D/YEST^ 0.1, analyze as drains on a barrier.
b. D/YEST - 0.8, analyze as drains above a barrier.
c. 0.1< D/YEST < 0.8, analyze both ways and use the
results with discretion. (The program analyzes this case as
drains above a barrier.)
Here, D is the distance from the horizontal centerline of the drain
to the barrier, and YEST is the maximum water table rise above the
drains at the midpoint between drains.
Although the criterion of (c) is not particularly satisfying, it is
a practical recommendation since the drainage system may mathematically
approach both limiting cases at different times. Improved accuracy
dictates use of the nonlinear solutions of Moody.(11)
The value of YEST is an input. Since the true value is unknown until
the analysis is completed, the input value is only an estimate (or
the design value when in the design mode). If an estimate is not
made and YEST is less than or equal to 0.01 foot (essentially zero),
the ratio D/YEST is presumed unknown. The value of D is then used to
select the case - if D is equal to or greater than 1.0 foot, analysis
is for drains above a barrier.
The program criteria can be overridden by means of the input variable
NTYPE. A value of 1 selects the case of drains above a barrier while
a value of 2 selects the case of drains on a barrier. Use of this
variable permits analysis by both cases in the intermediate range,
though this must be done by separate runs.
The flow chart of Figure 18 summarizes the selection process.
Note that it is possible to specify an analysis by the wrong method
according to either criteria (a) or (b). This situation is noted on
output, as is a recommendation that the intermediate range be analyzed
as both cases. The method of analysis, how it was selected, and the
recommended method of analysis based on the computed results are
printed on output, both for design and analysis.
120
-------
Anyother
DRAINS
ON
THE
BARRIER
FIGURE 18. PROGRAM CRITERIA TO SELECT DRAINAGE CASE
121
-------
APPLICATION OF METHOD
Two problems arise in using the USBR drainage methods:( )
a. Use of a linear solution for long drainout periods over-
estimates the water table decline.
b- Reinitializing the water table when an instantaneous
buildup occurs introduces volumetric errors because of the
mathematics involved.
The latter problem is insignificant for semiarid areas where deep per-
colation occurs at intervals of 10 days or more. It is serious for
humid areas where deep percolation from precipitation occurs fre-
quently, or when unsaturated flow results are used on a daily basis
(as when magnetic tape input via unit 1 is used). Errors in the mass
balance of the saturated system are the best indicators of this dif-
ficulty. The solution is to increase the short drainout periods by
combining several days of deep percolation into a single application.
The program implements this approach by means of the input variable,
YBLIM. If the deep percolation on a given day is not sufficient to
attain this value, it is accumulated and drainout proceeds. The next
buildup is added to this accumulated value and the sum is again com-
pared to YBLIM. The process continues and is repeated as necessary.
By using a larger value of YBLIM, the short drainout periods are
lengthened and mass balance usually improves. Rather than specifying
a minimum drainout period, the above approach was used because it is
more responsive to variations in the deep percolation inputs during a
run.
Since deep percolation inputs may produce either an instantaneous
buildup or decline, it is the absolute value of the accumulated deep
percolation that is compared to YBLIM. If YBLIM is left blank or
set to 0, the program uses an internally set value of 0.1 foot.
The first problem of long drainout periods is corrected by subdivid-
ing these periods into two portions. A zero buildup is assumed, the
water table is reinitialized, and the drainout process begins again.
There are two means available in the program to accomplish this:
a. The maximum allowable length of a drainout period in days can
be specified on input by MXTDR. If a drainout period reaches this
length, reinitialization occurs. If MXTDR is left blank or set to
0, the program uses a preset value of 60 days.
122
-------
b. A fictitious, near zero, deep percolation can be added on
the day when reinitialization is desired. This approach is use-
ful when input is by cards, but cannot be used with tape input.
The magnitude of the application must be sufficient to exceed
the value of YBLIM which, as noted above, cannot be set to zero
exactly. In general, if SY is the specific yield, the required
value is YBLIM x SY x 12 inches. When the value of YBLIM is
not important (the usual case when input is by cards), the val-
ues of YBLIM and MXTDR can be chosen to insure the desired result.
For example, if YBLIM = 0.0001 foot and SY = 0.1, a fictitious
application of (0.0001) (.1) (12) = 0.00012 inches is used on the
desired dates. The value of MXTDR is then set to a large value
(say 365) to assure it does not shorten the desired drainout
period.
DYNAMIC EQUILIBRIUM
If the same annual set of deep percolation inputs is used for succes-
sive years, the drainage cycle will approach a repeating annual pat-
tern. When this repeating cycle is attained, "dynamic equilibrium"
is said to have been achieved. It is this condition for which the USER
designs drains.
A convenient test for dynamic equilibrium is based on the comparison
of the water table position for successive years at the same time
within the annual cycle. Year-end values are practical for this pur-
pose, and the variable YTEST is read for the test. Its value can be
picked to obtain as strict a tolerance as desired. The smaller the
value, the more cycles required to satisfy the test.
In the program, testing for dynamic equilibrium does not occur until
all years of deep percolation inputs have been read (NYRS) and the
repetition mode is being used (NREPT > 0). Computations terminate
when the year-end water table positions are within YTEST.
DEEP PERCOLATION DATA
There is no restriction on the sign of the deep percolation inputs. A
positive value denotes water being added to the saturated system while
a negative one denotes water being removed. The latter case results
when unsaturated flow phenomena are explicitly used and evapotranspira-
tion produces upward flow.
The program allows deep percolation inputs on a daily basis. Deep per-
colation occurring on a given day is assumed to be applied at the end
of the day (at the start of the next day). For example, an application
on June 6 is applied as an instantaneous buildup or decline at the
start of June 7.
123
-------
Computations are performed day by day for an entire year. Partial
years are not allowed. Deep percolations read from cards include the
date on which they occur. All other days are assumed to have zero
inputs.
Tape input using unit 1 is on a daily basis, starting with some initial
day and ending with some final day. Days preceding the initial or fol-
lowing the final are assumed to have zero inputs. Organization of the
tape is shown in Figure 19. Because this tape is also used for plot-
ting unsaturated flow results, it contains information not needed by
the drainout program. In particular, only the Julian day and the cum-
ulative leachate are used.
The cumulative leachate refers to the accumulated amount of water that
was leached from the unsaturated into the saturated system. The drain-
out program uses the cumulative amounts on successive days to compute
the deep percolation by taking the difference. Since the drainout pro-
gram can be restarted and may be run using data starting with a year
other than the first, the initial cumulative leachate for computing
the first day's input is required. This is inputed by the value of
DCSTRT. Once the program is running, the cumulative values are car-
ried automatically.
The correct starting year is selected by means of the input variable
INFILE, described in the next section.
NEGATIVE STARTING POSITION
The water table position at the start of the initial year of the run,
YSTRT, may be specified as negative. A horizontal water table is
assumed as long as it is at or below the drain level.
Logically, use of a negative starting position is limited to the
analysis mode. The gradual buildup of the water table to the drain
level can then be followed. For design, a negative YSTRT merely
wastes computer time since the dynamic equilibrium condition is
required for design.
The following minor modifications to the computations and output are
noted:
a. Although mass balance results are printed, the water in stor-
age above the drains will be negative. It does not represent a
true quantity of water since the drain acts as the reference
level. However, the change in storage during the year is correct.
124
-------
B0T (beginning of tape,
load point)
-E0F (end of
file mark)
L
h i ghday
1 record/day, starting
ow day, ending at hi ghday
1 f i le/year , 1st year
2nd year
3rd year etc
Each dai Iy record i s
WRITE (2) I YEAR,
SF(J), U(J),
Z(J),
where
1 YEAR
I I
XT
CL
DEFAMC
Cl
HED
ETS
N
J
TN(J)
Z(J)
SF(J)
U(J)
s
i s
is
is
is
statement:
ETS, N, (TN(J),
written by the binary write
II, XT, CL, DEFAMC, Cl, HED,
J = 1,N)
the year number
the JuI ian day
the decimal part of the day (will always be
iO for end of day)
the cumulative leachate from unsaturated zone (cm)
the cumulative evapotranspirat ion that could
not be supplied to the crop because of low
moisture contents (cm)
is the cumulative infiltration at the surface (cm)
is the amount of water at the surface, remaining
to infiltrate (cm)
is the net amount of water which has passed
through the surface in either direction(cm)
is the number of nodes in the column
is the particular node number
is the volumetric mositure content (d imens i onless)
is the evapotranspiration not supplied since the
last output time on tape 5 (see unsaturated flow
program) (cm)
is the amount of water that has moved between
node J and J 1 since the last output time on
tape 5 (cm)
is the evapotranspirat ion rate (cm/day)
FIGURE 19. ORGANIZATION-DAILY MOISTURE TAPE
125
-------
b. The adjustment factor for daily discharges does not apply
when there is no discharge during the year. Therefore, the
message ADJUSTMENT FACTOR NOT REQUIRED appears.
c. Analysis of the drainage case is meaningless when the water
table is always below the drains. Therefore, both analysis and
output are deleted.
RESTART OF PROBLEM
A design problem cannot be restarted by this method. This section only
applies to analysis.
CLSTRT has been defined as the amount of deep percolation to be added
on the first day of the initial year of the run. As noted in the
previous section, deep percolation occurring on a given day is applied
at the start of the next. Therefore, an input on December 31 of the
previous year is carried over and applied on January 1. To allow for
this carryover on the initial year of the run, its value must be read.
The carryover for subsequent years is automatic. On completion of a
run, the value of CLSTRT for the next year is printed. It includes
the actual deep percolation of December 31 plus any accumulated deep
percolation not yet applied due to use of YBLIM. (See APPLICATION OF
METHOD, p. 122.)
The water table position at the end of the run is also printed. It
should be used as the value of YSTRT if the problem is restarted for
the next year.
When a problem is restarted, it is necessary to position the input and
output tapes on units 1 and 3 for the first operation if they are used.
This is accomplished by means of the inputs INFILE and INFIL3. INFILE
is the initial file to be read on unit 1 while INFIL3 designates the
first file to be written on unit 3. Since each file corresponds to a
single year's data, the word file is synonymous with year.
The program skips over (INFILE-1) or (INFIL3-1) files and positions
the tape at the start (first record) of the initial file. For the
first run of a problem, the file numbers may be set to 1 or left blank.
For subsequent runs, they must be set correctly. As an example:
a. Input tape 1 contains 3 years of data, with the third year
representing an average annual input to be used for subsequent
years.
b. The initial run computed results for 5 years.
126
-------
c. To restart the problem at the 6th year, INFILE = 3 and
INFIL3 = 6.
DESIGN METHOD
The determination of the proper drain spacing is a trial-and-error
process. An initial spacing is selected and the response of the
drainage system to the deep percolation inputs is analyzed. The
maximum water table rise at dynamic equilibrium is then compared
with the design value. If the computed rise is greater, the drains
must be spaced closer while a smaller computed rise means the drains
may be spaced farther apart.
Special considerations for design include:
a. Selection of the initial spacing.
b. Criteria for deciding when an adequate design has been found.
c. Method of adjusting the trial spacings to arrive at the
design spacing.
Each of these items is now discussed.
Design Method - Selection of Initial Spacing
There are two alternatives available to the program user:
a. The initial spacing can be specified on card 1 of group 2 as
SPACE. This may be based on experience, a rational analysis, or
pure guess.
b. The steady state Donnan formula can be used. This formula
makes use of a steady deep percolation rate (seepage across the
water table).
In practice, the steady-state seepage rate, SEEP, is usually based
on the deep percolation from one irrigation averaged over the time
between irrigations during the peak of the irrigation season. As
such, it represents a maximum rate at which water must be drained
in order to provide a sufficient drainage zone. Experience has
shown that when drains are on the barrier, this rate may be halved.
The steady seepage rate can be specified by the input variable SEEP
on the group 6 card. Alternately, the user may specify the starting
(IM0NST, IDATST) and ending (IM0NFI, IDATFI) inclusive dates for a
period over which the average seepage rate is computed. The total
127
-------
deep percolation for the period is used with the length of the
period to compute SEEP. If the starting date is not specified,
it is taken as January 1 while December 31 is used when the ending
date is not specified.
In summary:
a. If SPACE > 0, it is used as the initial spacing even though a
value of SEEP or the period dates are specified.
b. If SPACE - 0, the Donnan formula is used.
When the Donnan formula is used:
a. If SEEP * °» i1: is used-
b. If SEEP is - 0, the seepage rate is based on deep percolations
during the period IM0NST, IDATST through IM0NFI, IDATFI. If the
dates are not specified, the period is taken as January 1 through
December 31.
Design Method - Criteria of AdequateDesign
The design maximum rise is given as YEST on card 1 of group 2. If
the computed maximum rise at dynamic equilibrium is YMAX, the design
is deemed adequate when the absolute value of (YMAX - YEST) is less
than or equal to YCHECK. Here YCHECK is the tolerance read from the
group 6 card.
If YCHECK - 0, the program uses a default tolerance of 0.5 percent
of YEST. Thus, if YEST =4.0 feet, YCHECK =0.02 feet. Should
YCHECK < YTEST (the value used to determine when dynamic equilibrium
is reached), it is probable that an adequate design will not be found.
This is particularly true if YCHECK is much smaller than YTEST. To
avoid this possibility, the program resets YCHECK to YTEST if
YCHECK < YTEST.
The maximum number of trial spacings is specified as NTRYMX on the
group 6 card. If a satisfactory design is not obtained after NTRYMX
trials, a message is printed and execution terminates. When NTRYMX is
not specified, a program default value of 20 is used. This is also
the maximum value. Should the user specify NTRYMX > 20, the program
resets it to 20.
Design Method - Spacing Algorithm
When the first trial design is finished, the program uses the YCHECK
criteria to determine whether an adequate design has been found. If
128
-------
it has not, the spacing is either doubled or halved depending on the
results. If the spacing is doubled, it must also meet the requirement
that the new spacing be at least 100 feet. This procedure is contin-
ued until spacings both greater and smaller than the correct design
are found. With the correct spacing bracketed, the computed maximum
rises are used in conjunction with the design value and linear inter-
polation to obtain future trial spacings.
Normally, three trials result in the correct spacing. Thus, the
initial spacing is not critical.
INPUT CHECKS AND ERROR MESSAGES
A number of logical checks are made on the input data. Failure to
meet any of the tests results in an error message. These messages
are listed and explained in the next subsection. They are prefaced
by the notation:
******** ERRORS IN INPUT DATA ********
In addition, either of the following messages may appear when the
design mode is used:
* MAXIMUM NO. TRIALS REACHED, SPACING NOT FOUND
* DYNAMIC EQUILIBRIUM NOT REACHED, TERMINATE **
The first message results when the limit NTRYMX is reached without the
spacing being determined according to the YCHECK criteria. The second
message is obtained when dynamic equilibrium is not reached using the
YTEST criteria by NREPT years.
Should any of these errors be detected, the termination message:
******** EXECUTION TERMINATED ********
is written. When a normal termination is obtained after all computa-
tions are finished, the message:
******** COMPUTATIONS FINISHED. NORMAL EXIT. ********
is written.
OUTPUT
Program output consists of printed listings, punched cards, and com-
puted results written on magnetic tape.
129
-------
Printed data listed by the program on unit 6 is standard output.
It includes:
a. Input data read from cards including deep percolation inputs
when ICRDIN = 1.
b. A mass balance summary for each year of the run.
c. The volume of water removed by the drains for each day along
with monthly and annual values as well as the percent each month
is of the annual total.
d. The water table position for each day, as well as the mini-
mum and maximum rises and their dates.
e. The restart parameters YSTRT and CLSTRT for future years.
f. The method of analysis (case) used as well as the recommended
case.
Punched card output is written on unit 7 by the program. This out-
put is optional, being selected by NPUN =1. It consists of a title
card identical to the group 1 card on input followed by pairs of
cards containing monthly drain discharge volumes for each year of
the run (or dynamic equilibrium year for design). The volume cards
have the form
IYEAR, ICRD, (QM0N(I), I = 1, 6)
punched using an I5,I2,6E12.5 format. Here
a. IYEAR is the year number.
b. ICRD is the card order within the year. ICRD = 1 contains
data for the first 6 months and ICRD = 2 contains data for the
second 6 months.
c. QM0N(I) are the individual monthly volumes.
This card output is intended for use as input to the drain effluent
prediction program (DE0200) of the irrigation return flow quality
model.
Computed results are written by the program on unit 3 as optional
output when NSAVE = 1. Each year's output is one file on unit 3.
130
-------
Four logical records are written in each file as shown in Figure 20
These include:
a. Title information, the date of the computer run, and the
year number of the computed results.
b. Daily water table position.
c. Daily discharge volumes.
d. Monthly discharge volumes and the annual value.
All data is written without format in binary form. The magnetic
tape is intended as input to the utility program UT06 which plots
computed drainout results on a cathode ray tube.
Logical Error Messages
The following is a complete list of error messages:
EMI - DRAIN SPACING = ? LE 0.0
EM2 - DESIGN MAX RISE = ? LE 0.0
EMS - NO. YEARS OF DEEP PERC. INPUTS = ? MUST BE 1 FOR DESIGN
EM4 - DEPTH TO BARRIER = ? NEGATIVE
EMS - EFFECTIVE RADIUS = ? LE 0.0
EM6 - PERMEABILITY = ? LE 0.0
EM7 - STORAGE COEFFICIENT = ? LE 0.0 OR GE 1.0
EMI applies only to analysis, EM2 and 3 only to design, while the
remainder apply to both design and analysis. Each is discussed below:
EMI - The drain spacing SPACE must be specified as a nonzero
positive quantity for analysis.
EM2 - Drains cannot be designed unless the design rise YEST is
specified as a nonzero positive quantity.
EMS - Only a single design year of deep percolation inputs, NYRS,
is permitted.
EM4 - The drains must rest on or be above the barrier material;
hence, D must be zero or positive.
EMS - The drains must have a finite effective radius, ER.
EM6 - The permeability, PERM, must be a nonzero positive quantity.
EM7 - The storage coefficient, SY, must exceed zero if a transient
solution is possible and must not equal unity by definition.
131
-------
B0T (BEGINNING OF TAPE)
E0F (END OF FILE MARK)
oc
LOGICAL RECORD
LOGICAL RECORD 2 LOGICAL RECORD 3 LOGICAL RECORD 1
h- O
•< «t
LU LU
O_
LU Q£
ce o
J I
J I
J L
AIDENK I )
NDATE
IYEAR
QMQN( I )
QYR
J
(AIDENT (I), 1=1,10), (Y(l), 1=1,365) (Q(I),1=1,365 (QMQN(I), 1=1,12), QYR
NDATE, I YEAR
NOTATION
80 character alphanumeric title
date of the computer run
year number of computed results
water table position for Julian day I
daily drain discharge for Julian day I
monthly drain discharge for month I
annual drain discharge for year I year
FIGURE 2O. ORGANIZATION OF MAGNETIC TAPE OUTPUT
DRAINOUT SUBMODEL
-------
SECTION XII
USER'S MANUAL
INTERFACE FOR CHEMISTRY PROGRAM
INTRODUCTION
This program converts data generated by the Unsaturated Flow Sub-
model to forms suitable for input to the Unsaturated Chemistry and
Drainout Submodels.
The conversion for Unsaturated Chemistry involves a transfer from
a nodal concept in Unsaturated Flow to a segment concept in Unsatu-
rated Chemistry. Soil water flux and content values at each node
are transformed to (1) mean or representative flux values at the
upper and lower boundaries of each soil segment and (2) mean or
representative soil water contents for each segment (segment vol-
umes) . This procedure is repeated for each time step being used in
Unsaturated Chemistry. A commonly used time increment in 0.1 day.
At the end of each day, Interface can write binary information to
unit Tape 2. This information for Drainout is similar to the data
read from input unit Tape 3 with the exceptions that Q has been
added, and DEFAMC and K have been deleted. Additional information
can be found in the DATA DECK STRUCTURE and TAPE OUTPUT DETAILS
sections.
The current version of the Interface Program is designed to convert
data from 21 moisture nodes to data suitable for 10 soil segments
plus one segment to account for water applications at the surface.
A more generalized program could be developed at a future date.
However, it is possible to modify the program for other nodal and
segment combinations provided a few relationships are considered.
In general, a workable procedure is to determine the number of
soil segments needed based on dispersive considerations and then
compute the number of nodes to be used in Unsaturated Flow from
the equation:
No. Nodes = No. segments * N +1 [47]
Where:
N = 1, or 2 or 3 . . .
No. segments = the number of soil segments not counting
the segment for surface applications
133
-------
It is not recommended to use a value for N less than 1. Increasing
the magnitude of N increases the number of nodes, thus increasing
the accuracy and computer execution time of the Unsaturated Flow Sub-
model. Using 10 soil segments, values for N of 2 or 3 have proven
to be adequate.
The number of terms in the equation for segment volume given in the
Interface Program listing (card label, INTG 250) is given by
No. Terms = N +1 [48]
The first and last terms carry half the weight of the remaining ones.
The sum of the weighting fractions must equal the segment length in
cm. Segment length is computed by
Seg. Length (cm) = Length of Soil Column (cm) [49]
No. segments
The "DO Loop" parameters (Card label, INTG 230) must be set as follows:
First Parameter - equals 1
Second Parameter - equals N*(No. segments - 1) +1
Third Parameter - equals N
The "DO Loop" parameters (card label, INTG 280) must be set as follows:
First Parameter - equals N +1
Second Parameter - equals N*(No. segments - 1) +1
Third Parameter - equals N
Implied "DO Loop" termination in write (1) statement must be changed
to correspond to number of soil segments +1.
The card labeled INTG 350 must be set as follows:
MOISIN(M) = SF(M1)
where: M equals the number of segments +2.
Ml equals the number of nodes -1.
UNIT NUMBERS
Five file units are used by the Interface Program. Tape 3 is under-
stood to contain the binary output data written by the Unsaturated
Flow Submodel. Unsaturated Flow writes this data set out to Tape 5,
and Interface accepts it as input on Tape 3.
134
-------
Interface writes binary output to Tape 1 which is read by Unsaturated
Chemistry on unit Tape 1. Binary output data may be written to Tape 2
in a form suitable for input to the Drainout Submodel via its Tape 1.
In addition, use is made of the common input and output file units.
DATA DECK STRUCTURE
Data from two basic card groups are read from the common input file.
These groups contain title and control information.
CARD GROUP 1 - Title Information
Column
No.
Var Var
type name
Card 1 of 1
1-80 A
AID
Variable
description
Title or heading information used
to identify problem. Title is
printed on output.
CARD GROUP 2 - Control Information
1-5
6-10
IPRINT
11-15
16-20
21-25
26-30
31-35
36-40
41-45
I
I
I
I
I
I
I
INFILE (1
INFILE (2
INFILE (3
IRECB (1)
IRECB (2)
IRECB (3)
LSFILE
The number of moisture nodes used
in Unsaturated Flow Submodel.
Print interval for output listing
of data being written to Tape 1.
Value does not affect write inter-
val to Tape 1.
Initial file to use on Tape 1
Initital file to use on Tape 2
Initial file on unit Tape 3
First record initial file on
unit Tape 1
First record on initial file on
unit Tape 2
First record on initial file on
unit Tape 3.
Last file for unit Tape 3
135
-------
Column
No.
Var Var
typ_e_ name
46-50
51-55
56-60
IRECB
NOREC
NSAVE2
Variable
description
Last record for unit Tape 3
No. of records in full file on
unit Tape 3.
If = 1, Tape 2 is written in form
suitable for input to
Drainout Submodel.
If = 1, Tape 2 is not written.
TAPE OUTPUT DETAILS
Tape 1 - Output to Unsaturated Chemistry Submodel
WRITE (1) (IYEAR), II, INT, JDUM, SEGVOL (J) MOISIN (J), MOISONT
(J), TEN, U(J), J=l, 11)
Variable
name
IYEAR
II
INT
JDUM
SEGVOL (J)
MOISIN(J)
MOISOUT(J)
TEN
Variable
description
The year number.
The Julian day number.
The time step within each day.
A dummy variable reserved for
future use.
The volume of water (cm3/cm2) con-
tained in each soil segment, J.
Note: Segment SEGVOL (1) contains
water ponded on soil surface.
The flux of water (cm3/cm2) across
the upper boundary of segment J.
The flux of water (cm3/cm2) across
the lower boundary of segment J.
A dummy value for soil water tension
set equal to 300 cm water. Values
for TEN must be computed for each
segment in some applications of
Unsaturated Chemistry.
136
-------
Variable Variable
name description
U(J) The consumptive use of water
(cm3/cm2/day) for segment J.
Note: If uptake of NH4+ and
N03~ is taken to be propor-
tional to consumptive use,
then U(J) must be expressed
in terms of cm3/segment/time
step.
Tape 2 - Output to Drainout Submodel
WRITE (2) IYEAR, II, XT, CL, CI, HED, ETS, Q, (TN(J), Z(J), SF(J),
U(J), J=1,Q)
IYEAR The year number.
II The Julian day number.
XT The fraction of a day when the data
were written (e.g., 0.1 day, 0.2 day,
etc.)
CL Cumulative deep percolation (cm3/cm2)
CI Net infiltration from applied water and
errors in handling the upper boundary
(cm3/cm2).
HED Water standing at surface (cm).
ETS Total quantity of water crossing
surface boundary (cm3/cm2)
Q Number of moisture nodes.
TN(J) Volumetric moisture content (cm3/cm2)
at time XT for Node J.
Z(j) Deficit moisture (cm3/cm2) at Node J
during the last time interval written
to Tape 5 by Unsaturated Flow.
137
-------
Variable Variable
name description
SF(J) Total soil water flux (cm3/cm2)
between Node J and (J+l) during
the last time step written to
Tape 5 by Unsaturated Flow.
U(J) Evapotranspiration rate at Node J
(cm3/cm2/day).
138
-------
SECTION XIII
USER'S MANUAL
UNSATURATED CHEMISTRY PROGRAM
INTRODUCTION
The purpose of this manual is to provide the user with detailed oper-
ating instructions for the Unsaturated Chemistry Program. The manual
is written with the assumption that the user has read the sections of
this report and related references pertaining to the program's con-
cepts and theoretical basis. The manual includes specific discus-
sions of the types of data needed to run the program. Use or gener-
ation of all input and output media such as printed output, data cards,
magnetic tapes, and disks is discussed in detail.
The source program is written in FORTRAN IV computer language for a
Control Data Corporation CYBER 74-28 machine. The program as compiled
under the Fortran Extended Compiler requires 105K octal to load and
75K to execute.
CARD GROUPS
The basic card groups needed to run the program are illustrated in
Figure 21. Those groups marked with a dot must be omitted from the
deck when the nitrogen portion of the model is bypassed (NBYPAS = 1).
Group 6.5 must be included in the deck when the salt portion of the
model is bypassed (SBYPAS =1). At all other times, group 6.5 is
omitted.
In a similar manner, group 7 is included in the deck only when flow
data are to be read from cards (ITEST = 1). When ITEST = 0, flow
data are read from tape 1.
Group 8 contains the initial soil analysis and must be present when
the program is started from time zero (IRERUN = 0) . When IRERUN = 1,
group 8 is deleted from the deck.
Group 9, the restart data, is punched by the program at the conclu-
sion of a run provided IPUNCH = 1. Group 9 is read into the program
and must be present when IRERUN = 1 and IREADP =1. If IRERUN = 1,
but IREADP = 0, group 9 is not included, but tape 3 containing the
restart data must be present. This tape is written whenever
IPUNCH * 0.
139
-------
Group *I6 Crop
uptake of nitrogen
/Group* 15 Root
distribution
Group* 14 Organic
nitrogen applications
Group* 13 Fertilizer applications
,/Group *I2 Irrigation water
application dates
Group*ll Organic nitrogen
application dates
/Group* 10 Fertilizer
application dates
Group *9 Restart data
Group *8 Initial soil analysis
Group*? Flow data
/Group*6.5 Salt by-pass data
/ Group *6 Irrigation water .
analysis
Group *5 Temperature data
//Group*4 Constituent horizon
depths
Group *3 Temperature horizon
depths
Group *2 Control cards
/Group* I Title card
FIGURE 21. CARD GROUPS IN UNSATURATED
CHEMISTRY SUBMODEL
140
-------
CARD GROUP 1 - Title Card
Column
No.
Var Var
type name
Card 1 of 1
1-80 A TITLE
CARD GROUP 2 - Control Cards
Card 1 of 3
1-5 R DELX
6-10
21-25
26-30
31-35
36-40
41-45
R
R
R
X
R
DELT
11-15 X
16-20 R CH
CHI
AL
A2
PK
Variable
description
Any title desired on printed output.
Soil segment size (cm). See page 145 for
details.
Time interval size (days) corresponding
to intervals for moisture flow on tape 1,
Not used by program.
Convergence criterion for nitrogen trans-
formation subroutine. See page 145 for
details. Suggested values are in the
range 0.1 - 10.0 p/m nitrogen.
Convergence criterion for ion exchange
subroutine. See page 145 for details.
Suggested values are in the range
1 x 10"2 - 1 x 10~6 moles/liter.
Shutoff criterion for nitrogen trans-
formation subroutine. See page 145 for
details. Suggested values are in the
range 0.1 - 2.0 ug/soil segment/time
step.
Shutoff criterion for ion exchange sub-
routine. See page 146 for details.
Suggested values are in the range
1.0 - 5.0 ug/soil segment/time step.
Not used by program.
Fraction of total plant uptake of nitro-
•gen as N03-N. A suggested value is
0.95. As NOi-N mass -*• 0.0 plant uptake
of NOi-N -»• 0.0 regardless of plant
uptake imposed on system.
141
-------
Column
No.
46-50
Var Var
type name
R
PK1
51-55
R FACT
Card 2 of 3
1-5 I LL
6-10
11-15
16-20
21-25
26-30
31-35
36-40
MM
0
CR0P
T0
NT
IPRINT
JPRINT
Variable
description
Fraction of total plant uptake of nitro-
gen as NH^-N. A suggested value is
0.05. Fraction as N03-N plus fraction
as NH^-H must equal unity. The same
relationship for N03-N holds for NHH-N
uptake with respect to the zero mass
boundary condition.
Plant-N uptake constant for use with
consumptive use values (if used).
Starting day of run relative to time
reference date (e.g., if reference date
is January 15 and starting day for run
is January 21, the relative starting
day would be day 7).
Termination day of run relative to time
reference date.
Number of soil chemistry horizons (must
be in range 1 to 9 inclusive).
If 1, consumptive use used to determine
N-uptake by plants.
If 2, reserved for future use.
If 3, read plant N-uptake data from
cards (see card groups 15 and 16).
Number of temperature horizons (must be
in range 1 to 4 inclusive).
Number of weekly temperature data cards
(must = 0 for NBYPAS = 1) (see group 5),
Print interval (days) for predicted soil
chemistry data.
Time interval within a day on which the
IPRINT and IMASS print options are to
be activated.
142
-------
Column
No.
41-45
46-50
51-55
56-60
61-65
66-70
71-75
Var Var
type name
I INK
I IRERUN
I IPUNCH
76-80 I
Card 3 of 3
1-5 I
6-10 I
IREADP
ITEST
START
SM0NTH
YEAR
ISTOP
IMASS
Variable
description
Write interval (days) for output to
tape 2 (see section 3).
If 0, program is to be run from initial
input deck.
If 1, program is to be restarted from
tape or restart deck.
If 0, no restart data deck will be
punched (a restart tape (No. 3) is
written).
If 1, a restart data deck will be
punched.
If 0, program is to be restarted from
tape 3 or started from initial data.
If 1, program is to be restarted from
cards (see group 9).
If 0, moisture flow data are read from
tape 1.
If 1, moisture flow data are read from
cards (see group 7).
Reference date in reference month (e.g.,
March 21 would be STARTING DAY =21).
See below.
Reference month. Month to which all
starts of the program are referenced
(e.g., if reference month is March,
STARTING MONTH = 3).
Year for this run (e.g., if this is
the third year being run, YEAR = 3).
Year when run terminates (e.g., if run
is to terminate after 3 years and starts
with year 3, ISTOP =5).
Print interval (days) for summary of
nitrogen balance for system.
143
-------
Column
No.
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
Var Var
type name
I IPRINTI
I IPRINTJ
I NBYPAS
I SBYPAS
I IDYSTR
I IDYSTP
I INFIL2
I INFIL1
I ICONT1
I ICONT2
I IOP
Variable
description
If 1, basic control card data are
printed.
If 0, printing of above is suppressed.
If 1, input data are printed.
If 0, printing of above is suppressed.
If 0, nitrogen transformations are
considered.
If 1, nitrogen transformations are
bypassed, salt reactions remain in
model if SBYPAS = 0.
If 0, salt reactions are considered.
If 1, salt reactions are bypassed, nitro-
gen transformations remain in model if
NBYPAS = 0.
Starting day for run. If LL > IDYSTR,
LL - IDYSTR Records are skipped on
tapes.
Stopping day for run (uses MM value if
year = ISTOP).
Starting file on output tape 2 (used
when making a restart on tape 2 con-
taining previous files).
Starting file on input tape 1 (used
when more than one year's data are to
be inputed to flow calculations).
If
If -
If =
If
1, a new set of water application
dates is read in each year.
Also, no records can be skipped
on tape 1.
1, no record skipping can occur on
tape 2. Otherwise (LL - IDYSTR)
records are skipped on tape 2
following file positioning.
1, lime precipitation is checked
in relation to available pore
space.
0, option is bypassed.
144
-------
Column
No.
65-70
IOPN
Variable
description
If = 1, transition-state nitrification
subroutine used.
If = 0, statistical nitrification
equation.
(If = 1, constants and one statement
function must be inserted
into subroutine NITR; details
are listed in the subroutine
as comments.)
71-75
IPC02
75-80
IREK
If
If
If =
If =
ADDITIONAL DETAILS FOR CARD GROUP 2
1, user supplies pC02 valves used
(see group 8).
0, program relates pC02 to soil
moisture content.
1, cation exchange coefficients
for Ca-Na and Ca-Mg ion
exchange are read from soils
analysis data cards.
0, default values within program
are used.
Soil Segment Size
Soil segment size (or DELX, see figure 1) must be in the range [depth
(cm) to water table/(25 - 1)] <_ Soil Segment Size <^ the size (cm) of
the largest chemistry or temperature horizon. Here the 25 refers to
array size currently used in the program.
Converge1
This constant determines the accuracy of the nitrogen changes computed
in subroutine TRNSFM. A smaller value will increase accuracy and pro-
gram execution time.
Converge2
This constant determines the accuracy of mass changes in the system
computed in subroutine XCHANGE. A smaller value will increase accuracy
and program execution time.
Check1
This constant is used in subroutine CHK to determine if subroutine
TRNSFM should be called for that time step or bypassed and the
145
-------
previously computed nitrogen changes used. If the rates of change
for nitrogen are less than CHECK1, TRNSFM is bypassed. However,
TRNSFM is called at least once each day to avoid possible "drift"
in the calculations.
Check2
This constant is used in subroutine CHK to determine if subroutine
XCHANGE should be called for that time step or bypassed, in which
case previously computed changes due to ion exchange, etc., are set
equal to zero. XCHANGE is called at least once each day.
CARD GROUP 3 - Temperature Horizon Depths (Omit if NBYPAS = 1)
Column Var Var Variable
No. type name description
Card 1 of 1
1-5 R TH0R(J) Depth (cm) from surface to bottom of
first temperature horizon.
6-10* R TH0R(J) Depth (cm) from surface to bottom of
second temperature horizon (if it
exists).
(J = subscript for each horizon.)
CARD GROUP 4 - Chemistry Horizon Depths
Card 1 of 1
1-5 R H0R(J) Depth (cm) from surface to bottom of
first chemistry horizon.
6-10** R H0R(J) Depth (cm) from surface to bottom of
(J+l) second chemistry horizon.***
(J = subscript for each horizon.)
* The program allows the inclusion of up to four temperature horizons.
If included, the depths for numbers 3 and 4 would be punched in col-
umns 11-15 and 16-20, respectively.
** The program allows inclusion of up to nine chemistry horizons. If
included, the depths for numbers 3, 4, 5, 6, 7, 8, and 9 would be
punched in columns 11-15, 16-20, 21-25, 26-30, 31-35, 36-40, 41-45,
respectively.
*** If only one chemistry horizon is desired, only one need be included.
146
-------
CARD GROUP 5 - Temperature Data (Omit if NBYPAS = 11
Column Var Var Variable
No- tyP6 name description
Card 1 of NT*
*~ 2 Card indentification (not read by
program).
3-10 R TT(K) Weekly temperature (°C) for uppermost
(first) temperature horizon.
10-20** R TT(K) Weekly temperature (°C) for second
(K+l) temperature horizon.***
(K = subscript for weekly temperatures.)
CARD GROUP 6 - Irrigation Water Analysis
Card 1 of 1
(All units in meq/L of the species
shown.)
1- 5
6-10
11-15
16-20
21-25
26-30
R
R
R
R
R
R
ANH3(1)
AN03(1)
CA(1)
ANA(l)
AMG(l)
HC03(1)
MV
N03"
Ca++
Na+
Mg++
HC03
* NT is the number of temperature data cards. Each card contains
weekly data for all temperature horizons.
** The program allows inclusion of up to four temperature horizons.
If included, temperature data for horizons 3 and 4 would be punched
in columns 21-30 and 31-40, respectively.
*** If only one temperature horizon is desired, only one need be
included.
Note: The first temperature data should correspond to the week begin-
ning" with the reference data (see Section 5 for exceptions).
147
-------
Column Var Var Variable
No. type name description
31-35 R CL(1) Cl"
36-40 R C03(l) C0=3
41-45 R S04(l) 50%
Note: All irrigation water is assumed to have same analysis.
CARD GROUP 6.5 - (Needed if SBYPAS = 1, and NBYPAS = 0)
Card 1 of Q-l
1-5 R SERATI0(J) Fraction of total NHi,* which is soluble
NH,/.
6-10
11-15
CARD GROUP 7
Card 1 of NS
1-10
11-20
21-30
31-40
41-50
R
R
-
*
R
R
R
R
R
U(J)
ACTCA(J)
Moisture Flow
if ITEST = 0,
CMH201 (J)
M0ISIN(J)
M0IS0UT(J)
TEN(J)
u(J)
Ionic strength - needed only when
IOPN = 1.
Ca"1"* activity - needed only when
IOPN = 1.
Data (Needed only if ITEST = 1;
tape 1 is used)
Volume of water in segment (cc) .
Water movement across upper boundary
of segment (cc/time step).**
Water movement across lower boundary
of segment (cc/time step).**
Average moisture tension (cm) over
time step (tensions shown as negative)
Plant water uptake (cm/time step/soil
segment) .
* NS is the number of moisture flow data cards. There must be a card
for each segment including No. 1. For details in computing NS, see
page .
** Positive sign for downward flow, negative for upward flow.
148
-------
CARD GROUP 8 - Initial Soil Analysis (Omit if IRERUN = 1)
Column
No.
Card 1 of
Var
type
NC*
Var
name
Variable
description
(Units in meq/L of species shown unless
otherwise noted.)
1- 5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
66-70
71-75
76-80
R
R
R
R
R
R
R
R
R
R
R
R
R
R
R
R
ANH3(1)
AN03(1)
UREA(l)
CA(1)
ANA(l)
AMG(l)
HC03(1)
CL(1)
C03(l)
S04(l)
EC(1)
XX5(1)
CAL(l)
BD(1)
SAMT(l)
CNl(l)
NHi4 (soluble form) .
NO 3".
UREA.
Ca+*
Na+.
Mg*+.
HC03".
Cl".
C03=.
SO,,'.
Exchange capacity (meq/100 gm
Gypsum (meq/100 gm) soil.
soil) .
If 0.0, no lime is present in soil.
If 1.0, lime is present in soil.
Bulk density (gm/cc of soil) .
Soil organic nitrogen (ug/gm
Carbon-nitrogen ratio of soil
matter (e.g., 40:1 would be
as a 40.0).
soil) .
organic
punched
NC is the number of chemistry horizons.
149
-------
Column Var Var Variable
No. type name description
Card 2 of Set
1-5 R XTRCT(l) Gm water/gm soil in soil extract.
6-15 R PC02(1) C02 partial pressure in atmospheres
(needed only if IPC02 = 1)
16-20 R AKCS(l) Ca-Na exchange coefficient.
21-25 R AKCM(l) Ca-Mg exchange coefficient.
CARD GROUP 9 - Restart Data
needed only
if IREK = 1
This data group is punched by the program during a previous run. See
page for details. The punch/read format is (4I5/, (6E13.5/6E13.5/
6E13.5/6E13.5/E13.5)).
CARD GROUP 10 - Fertilizer Application Dates
Card 1 of 1*
1-5 I IT0T Number of fertilizer applications.
6-10 I IADD(K) First application day relative to ref-
erence starting date.
11-15** I IADD(K) Second application day relative to ref-
(K+l) erence starting date.
CARD GROUP 11 - Organic Nitrogen Application Dates (Omit if NBYPAS =1)
Card 1 of 1
1-5 I JT0T Number of organic nitrogen applications.
6-10 I I0RNAP(K) First organic nitrogen application day
relative to reference date.
* A second card may be used if needed (i.e., if the number of appli-
cations exceeds 15, in which case the 15 format still applies).
** The program allows inclusion of up to 25 fertilizer application
dates. The same format, 15, is used for any additional dates desired.
150
-------
Column Var Var Variable
No. type name description
H-15* I I0RNAP(K) Second organic nitrogen application day
(K+l) relative to reference date.
CARD GROUP 12 - Irrigation Water Application Dates
Card 1 of 1**
1-5 I IRT0T Number of irrigations.
6-10 I IRR(K) First irrigation day relative to ref-
erence date.
11-15*** I IRR(K) Second irrigation day relative to ref-
(K+l) erence date.
CARD GROUP 13 - Fertilizer Applications
Card 1 of NF****
All units in Ibs/acre of species shown
unless otherwise labeled.
1-5 R FERT(l) Depth (cm) of a uniform fertilizer
application.
If = 0.0, a surface application is
indicated.
6-10 R FERT(2)
11-15 R FERT(3) N03".
16-20 R FERT(4) Urea (NH2-C-NH2).
* The program allows inclusion of up to organic nitrogen application
dates. If included, application dates for applications 3, 4, and 5
would be punched in columns 16-20, 21-25, and 26-30, respectively.
** A second data card may be needed if more than 15 irrigation dates
are us ed.
*** The program allows inclusion of up to 25 irrigation dates. If
included, applications 3-25 would follow the same 15 format.
**** NF is the number of fertilizer applications. The maximum number
allowed in the program is 25.
151
-------
Column Var Var Variable
No . type name description
21-25 R FERT(5) Ca**.
26-30 R PERT (6)
31-35 R FERTC7) C03=.
CARD GROUP 14 - Organic Nitrogen Applications (Omit if NBYPAS = 1)
Card 1 of NO*
1-5 R 0FERT(1) Depth (cm) of uniform organic nitrogen
applications (surface applications are
not allowed) .
6-10 R 0FERT(2) Carbon: nitrogen ratio of organic matter
added (containing the organic N) (e.g.,
40:1 would be punched as a 40.0).
11-15 R 0FERT(3) Organic nitrogen added (Ibs/acre) as
oven dry weight.
CARD GROUP 15 - Plant Root Distribution
Card 1 of 1
1-10
11-20
21-30
31-40
41-50
51-60
R
R
R
R
R
R
KP2(1)
KP2(2)
KP2(3)
KP2(4)
KP2(5)
KP2(6)
Fraction (decimal) of
0-1 foot depth.
Fraction
Fraction
Fraction
Fraction
Fractions
from 1 foot
from 2 feet
from 3 feet
from 4 feet
> 5 feet.
plant roots from
to 2
to 3
to 4
to 5
feet.
feet.
feet.
feet.
* NO is the number of organic matter applications. The maximum num-
ber allowed in the program is 5.
Note: Plant root distribution is independent of time. If roots do
not extend to 5 feet, use zero fractions for roots below their deepest
extension.
152
-------
CARD GROUP 16 - Plant Uptake of Nitrogen
Column
No.
Var
type
Var
name
Variable
description
Card 1 of NP*
1-10 X
11-20 R
UPTK(K)
Space for card identification.
Plant uptake of nitrogen (lbs-N/acre/
semimonthly).
(K = subscript for each semimonthly
period.)
TAPE INPUTS
Tape input to the Biological-Chemical program consists of two optional
tapes written in binary. Tape 1 contains moisture flow data written
by program INTFACE. An alternate procedure is to input moisture data
on cards. In this case a single set of constant moisture data is
read. Tape 3 may be read to restart the program after a previous run.
This tape is always written at the conclusion of a run. Unless a mag-
netic tape is equipped to tape unit 3, the binary data are written on
a disk or drum. Disk or drum data may be lost at the conclusion of
a run. The alternative to a restart from tape 3 is a restart from
data cards.
TAPE 1
Var
type
Variable
description
Logical Record
1 of NR
1 I
II
Reserved for future use.
* NP is the number of semimonthly** plant uptake data cards. The
program requires that 24 cards be present, even if not all of them
are used.
** Semimonthly means the period from the 1st day through the 15th
day or from the 16th day through the last day of the month. The
1st data card begins with the semimonthly period containing the
reference starting date (see page 143). Remaining data cards con-
tain data for consecutive semimonthly periods taken chronologically
from the 1st.
153
-------
Field
No.
2
3
4
5
Var
type
I
I
I
R
Var
name
12
13
13
CMH201(J)
Reserved
Reserved
Reserved
Current
for
for
for
Variable
description
future use.
future use.
future use.
column of water
8 R
9 R
TAPE 3
Logical Record
1 of 1
1 I
2 I
M0ISIN(J)
M0IS0UT(J)
TEN(J)
U(J)
IC0UNT
INFERTIN
N0RGRIN
NTEMPIN
Jth* soil segment .
Moisture flow (cc/time step) into the
Jth soil segment.
Moisture flow (cc/time step) from the
Jth soil segment.
Average current moisture tension (cm
for Jth segment (sign - ±) .
Consumptive use (cc/time step) for the
Jth segment .
Restart counter for plant uptake data.
Restart counter for fertilizer applica-
tion data.
Restart counter for organic matter appli
cations data.
Restart counter for temperature data.
* This sequence of nine fields is repeated Q times per logical
record. Q is the total number of segments. Its value may be
computed using the algorithm outlined on page 155. NR is the
number time steps for which data are written on tape 1.
154
-------
Field
No.
5*
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Var
type
R
R
R
R
R
R
R
R
R
R
R
R
R
R
R
R
R
Var
name
ANH3 (J)
AN03 (J)
UREA(J)
CA(J)
ANA(J)
AMG(J)
HC03 (J)
CL(J)
C03(J)
S04(J)
EC(J)
XX5(J)
CAL(J)
BD(J)
SAMT(J)
CN1 (J)
ORN(J)
Variable
description
pg of NHij-N contained in Jth** segment
pg of NOa-N contained in Jth segment.
pg of Urea-N contained in Jth segment.
pg of Ca contained in Jth segment.
pg of Na+ contained in Jth segment.
pg of Mg++ contained in Jth segment.
pg of HC03" contained in Jth segment.
pg of Cl~ contained in Jth segment.
pg of C03~ contained in Jth segment.
«•
pg of SOit contained in Jth segment.
Exchange capacity of Jth segment.
Gypsum concentration in Jth segment
(moles/L) .
If 1.0, lime is present in Jth segment
If 0.0, lime is not present in Jth
segment .
Soil bulk density of Jth segment.
pg of organic matter just applied to
Jth segment.
C:N ratio of organic matter in Jth
segment .
pg of organic-N in Jth segment.
* Fields No. 5-33 inclusive are repeated Q times in each logical
record. Q is the number of segments and may be determined using the
following algorithm: Q = DEPTH to W.T./DELX + 1.1, where any frac-
tion is dropped in determining the final value for Q.
** The segments begin with J = 1 and end with J = Q.
155
-------
22
23
24
25
26
27
28
29
30
31
32
R
R
R
R
R
R
R
R
33
CARD OUTPUTS
RN(J)
RC(J)
E5(J)
C5(J)
SA5(J)
CAS0(J)
AGS0(J)
BNH4(J)
XTRCT(J)
ANETLIM(J)
AZE(J)
IIK(J)
Residue nitrogen in Jth segment (pg/g).
Residue carbon in Jth segment (pg/g).
Exchangeable Ca
(moles/g).
Exchangeable Mg
(moles/g).
in Jth segment
in Jth segment
Exchangeable N in Jth segment
(moles/g) .
Undissociated CaSOit in Jth segment
(moles/L) .
Undissociated MgSO^ in Jth segment
(moles/L) .
Exchangeable NH
(moles/g) .
in Jth segment
Ratio of water/ soil in the extract for
Jth segment .
Current amount of pore space occupied by
lime in Jth segment.
A value proportional to K1 and asso-
ciated with lime precipitation-
dissolution and CC>2 partial pressures.
Value is for Jth segment.
Current lime indicator for Jth segment.
Card output from the program consists of data needed to restart a
run after a previous execution. This card deck is generated only
if the proper control card request has been made (see Card Group 2).
No card deck is punched if a program run terminates abnormally (e.g.,
a fatal error occurs). A discussion of the punch format is presented
in Section 1, under card group 9. See tape 3, page 154, for variable
descriptions.
156
-------
TAPE OUTPUTS
Tape output from the Biological -Chemical program consists of one
binary (tape 2) and one binary (tape 3) file. Tape 2 contains the
soil leachate data. This may serve as input to routines which con-
sider saturated flow and saturated chemistry. Tape 3 contains restart
data which may be read into the program to restart a previously ter-
minated run.
TAPE 2
Logical
Record
1 of NR
1
2
3
4
5
6
7
8
9
10
11
12
13
14
Var
type
I
I
R
R
R
R
R
R
R
R
R
R
R
R
Var
name
YEAR
IDAY
SUM0UT
DELN
AMT(l)
DEL(l)
AMT(2)
DEL(2)
AMrr(3)
DEL(3)
AMT(4)
DEL (4)
AMT(5)
DEL(5)
Variable
description
Year number.
Day number.
Water leached since start of run (cc/cm2).
Water leached since previous record was
written (cc/cm2).
N03 N leached since start of run (yg/cm2).
N03 N leached since previous record was
written (yg/cm2).
NHi/N leached since start of run (yg/cm2).
NH^~N leached since previous record was
written (yg/cm2).
Urea-N leached since start of run (yg/cm2)
Urea-N leached since previous record was
written (yg/cm2).
Ca** leached since start of run (yg/cm2).
Ca*"1" leached since previous record was
written (yg/cm2).
Na* leached since start of run (yg/cm2).
Na* leached since previous record was
written (yg/cm2).
157
-------
Logical
Record
15
16
17
18
19
20
21
22
23
24
Var
type
R
R
R
R
R
R
R
R
R
R
Var
name
AMI (6)
DEL(6)
AMI (7)
DEL (7)
AMT(8)
DEL(8)
AMI (9)
DEL (9)
AMT(IO)
DEL(IO)
Variable
description
Mg leached since start of run (yg/cm2) .
Mg leached since previous record was
written (yg/cm2) .
HCO~3 leached since start of run (yg/cm2)
HCO 3 leached since previous record was
written (pg/cm2) .
Cl~ leached since start of run (yg/cm2).
Cl~ leached since previous record was
written (pg/cm2) .
C0~3 leached since start of run (pg/cm2).
C0~3 leached since previous record was
written (yg/cm2).
SO", leached since start of run (yg/cm2) .
S0~. leached since previous record was
written (yg/cm2).
NR is the number of records written on tape 2. The write interval
(i.e., time between write operations is specified by the user (see
card group 2)).
HINTS ON PROGRAM USE
a. Be sure core is set to zero before program is loaded.
b. The program can be operated very easily from a time-share ter-
minal. This should be done when possible.
c. Subroutines will not load in a random order due to varying sizes
of common blocks. See suggested load sequence on page 160.
d. An end-of-file is written on tape 2 at the conclusion of a run
and at the end of each year. This means that output from several
runs and/or years may be stored as separate files on the same or
individual tapes.
158
-------
e. A run which exits due to a time limit or other error will not
produce a restart deck or tape.
f. The number of temperature cards may be reduced when making a
restart run by altering the value of the number contained in
Field No. 4 of the restart deck or tape. The number placed there
by the program is the number of cards or records to skip in order
to locate the correct temperature data card.
g. If no plant uptake of N is required, set CROP equal to 1 and
leave columns 51-55 of the first control card blank. In this
case, CARD GROUP 16 may be omitted from the deck.
h. When using the IOPN = 1 option, values for certain parameters,
plus a statement function must be inserted into NITR.
i. The program can be operated with either the salt portion or the
nitrogen portion bypassed.
j . The program can be operated so that certain months of the year
are bypassed.
k. Be sure to check the input data carefully. The program may run
and terminate normally with bad data.
1. The mass balance check on nitrogen is an important item. Any
large errors here indicate something has gone wrong.
m. If a run should "blow up" first check your card input data and
then any data being read from tapes.
n. An average execution time on a CDC CYBER 74-28 is about 0.4
second/day using 11 segments and 0.1 day time steps.
159
-------
SUGGESTED LOAD SEQUENCE
Load order
Routine
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
MAIN
COMBINE
TRNSFM
XCHANGE
EQEXCH
NITR
NEGN
SALTBP
EXECUTE
MCHECK
OUTPT
DAY
I DAY
THEDATE
UNITS 1
FL
PRNT
WGT
UPTAKE
CHK
TEMP
QROTP
QTRAN
COMP
PHCALC
SKIP
160
-------
SECTION XIV
USER'S MANUAL
SATURATED FLOW PROGRAM
INTRODUCTION
This program was designed to predict the average or steady-state
response of a subsurface drainage system to deep percolation inputs.
The drainage system consists of parallel, equally spaced, tile drains.
Either the average annual deep percolation rate for uniform seepage
(boundary condition 1, BC1) or the average annual water table position
for the ponded water case (boundary condition 2, BC2) provide the
flow inputs. In particular, this program was developed as part
of an irrigation return flow water quality model.
For computation purposes, the saturated aquifer is divided into a
number of individual stream tubes. Discharge rates through each
tube are identical. In addition, a node concept is employed to
permit the determination of each tube's location. Theoretical
techniques and details are contained in the program description.
Output includes the position of the water table (BC2) or the seepage
distribution between drains (BC1), location of the stream tubes, time
of travel down each tube, pore volume of each tube, and the steady-
state discharge from the drains. Pore volumes are intended for input
to the saturated chemistry portion of the water quality model while
the drain discharge and travel times are input to the drain effluent
prediction portion.
UNIT NUMBERS
Two input-output devices are referenced by the program:
Unit Unit description
designation or purpose
5 The common input tape or unit (CIT). It is the ,
computer system's standard input device.
6 The common output tape or unit (C0T). It is the
computer system's standard output device.
In addition, the program uses an internal scratch or peripheral
storage unit (unit 8) which is assigned to disk. Data defining the
streamline locations are written without format (binary) if a cathode
161
-------
ray tube (CRT) plot of the streamlines using a distorted vertical
scale is requested.
DATA DECK STRUCTURE
Seven input cards are required for BC1 and eight for BC2 as shown
in Figure 22. Certain cards are common to both boundary conditions.
Referring to the card numbers of Figure 22, the input deck for BC1
consists of cards 1, 2, 3, 4A, 5A, 7, and 8 while for BC2 cards 1,
2, 3, 4, 5, 6, 7, and 8 are used.
CARD GROUP 1 - Title Information
Column
No.
Card 1 of 1
1-80
Var
type
A
Var
name
TITLE j
Variable
description
An 80-character alphanumeric title
used to identify all printed and
plotted output. If this title is
desired to be centered on output,
it must be centered on the card.
CARD GROUP 2 - Problem Selection
Card 1 of 1
1-40 I
(4011 Format)
IPR0B(I) A control array which calls the proper
overlay for execution. Individual
programs in the water quality model
were originally designed for inclusion
in an overlay system. Presently
(March, 1973), only the saturated flow
and drain effluent portions are
included. When IPR0B(I) is set to 1,
the Ith overlay is called. Thus
IPR0B(1) = 1, call the saturated
flow routine.
IPR0B(2) = 1, call the drain efflu-
ent routine.
On output, values of IRP0B(I) are
listed as I0C values under the
heading - PROBLEM CONTROL INFORMATION
FOR MAIN PROGRAM.
162
-------
No. of stream
tubes
/Grid division
- 1
/Geometry
t
\
. - .- , -A—
^~~
(4 A
^ -•
\
Boundary input/output
options
Problem selection
/TitU
Porous media
properties
FIGURE 22. DATA DECK STRUCTURE-SATURATED FLOW SUBMODEL
163
-------
Column
No.
Var
type
Var
name
Variable
description
CARD GROUP 5 - Boundary Condition and Input/Output Options
Card 1 of 1
1-10
IBC
11-50
(4011 Format)
Boundary condition selector.
IBC = 1, selects the ponded water
case.
IBC = 2, selects the uniform seepage
case.
IBC = any other value will produce
questionable results.
(There is no check for a
legitimate boundary con-
dition number.)
The boundary condition number is
incorporated into a subtitle on
output.
I0C(I) = An input/output option array
used to select printed and
plotted output. (See
Input/Output Options.)
CARD GROUP4 - Geometryfor BC2 (see Figure 25)
Card 1 of 2
1-10
R
DG
11-20
21-30
Card 2 of 2
1-10
R
R
R
A
ER
DG
11-20
AG
Distance from horizontal centerline
of drain to the impermeable barrier
(feet).
Drain spacing (feet).
Effective drain radius (feet).
Geometry for BC1 (see Figure 23).
Distance from horizontal centerline
of drain to the impermeable barrier
(feet).
Drain spacing (feet).
164
-------
-Water table
• Perm • / .'
. pore
o
o • .
- • . • • "-.-•..'• 17 .
-Impermeable barrier .' '. • • • • y
wr^//&ff^//j^^^^^ff^f^//^
AG
(a) BOUNDARY CONDITION 1
Uniform seepage
• ^ Impermeable barrier • .o
cO " ' . " ' ' 4 • . ' ,
Seep
(b) BOUNDARY CONDITION 2
FIGURE 23. BOUNDARY CONDITIONS-
SATURATED FLOW SUBMODEL
165
-------
Column Var Var Variable
No. type name description
21-30 R H Vertical distance from the water table
to the horizontal centerline of the
drain (feet).
31-40 R ER Effective drain radius (feet).
CARD GROUP 5 - Grid Division Data for BC2
Card 1 of 2
1-5 I Ml The number of primary grids to be
subdivided in the horizontal
direction.
6-10 I Nl The number of primary grids to be
subdivided in the vertical
direction.
11-15 I L The number of subdivisions of the
primary grid.
16-25 R HG The nominal size of the primary grid
(feet).
Card 2 of 2 Grid division data for BC1.
See above Ml Same description as in card 5.
Nl
L
CARD GROUP 6 - Seepage Rate for BC2
Card 1 of 1
1-10 R SEEP Average annual seepage or deep per-
colation rate between the drains
(feet3/feet2/year or feet/year).
CARD GROUP 7 - Porous Media Properties
Card 1 of 1
1-10 R PERM Coefficient of permeability (feet/year)
166
-------
Column
No.
11-20
Var
type
R
Var
name^
P0RE
Variable
description
Total porosity expressed as a decimal
(feet3/feet3 or dimensionless).
CARD GROUP 8 - Number of Stream Tubes
Card 1 of 1
1-10
STANDARD OUTPUT
NST
The number of stream tubes or stream-
lines used to define flow from one
side of the drain.
Output includes results printed on unit 6 and plotted on a cathode
ray tube (CRT) plotter*. Certain portions of the printed output
are standard while others are optional. All CRT plots are optional.
The next section describes the optional output.
Standard printed output common to both boundary conditions includes:
a. input/output options (I0C(I));
b. boundary condition number;
c. given and modified aquifer geometry;
d. porous media properties;
e. grid division data;
f. number of stream tubes;
g. flow into the drain per foot of length;
h. travel time down each tube;
i. pore volume contained in each tube;
j. error messages, if any.
In addition, standard output for BC1 includes the average seepage
rates between adjacent streamlines as well as the average rate
between drains. For BC2, additional output includes the ratio of
the seepage rate to the permeability and the maximum and minimum
water table position at 50, equally spaced points between the drain
and the midpoint between drains.
* Plotter subroutines are not included in the model version included
in this Volume.
167
-------
INPUT/OUTPUT OPTIONS
Individual output options are selected by setting the appropriate
I0C(I) value to 1. Options are bypassed if the value is left blank
or set to 0. Available options are identical for both boundary condi-
tions with the exception of I0C(28) which only applies to BC1. The
descriptions below may be more meaningful if reference is made to
the computational approach described in the program description. It
should be noted that most of the print options are intended for
debugging use.
I_ Option
1 Select all print options. All individual print options
(1=3 through 25) are set to 1 by this option.
2 Select all CRT plot options.* All individual plot options
(1=26 through 40) are set to 1 by this option.
3-25 Reserved for individual print options.
3 Print required values of streamlines used to define the
flow system.
4 Print coordinates of all nodes.
5 Print coordinates and the stream function at all nodes
that are required to define the system.
6 Print coordinates and corresponding potentials for all
points locating the required streamlines after horizontal
interpolation.
7 Print coordinates and corresponding stream function for
all nodes in each vertical column of nodes.
8 Print coordinates and corresponding potentials for all
points locating the required streamlines after horizontal
and vertical interpolation.
9 Print coordinates and corresponding potentials for all
points defining a streamline before sorting (raw data).
10 Print coordinates and corresponding potentials for all
points defining a streamline after sorting (sorted data).
* Not available in this version of the model.
168
-------
I Option
11 Print coordinates and corresponding potentials for all
points defining a streamline after the sorted data have
been modified to exclude points within the drain and a
point added on the drain perimeter (modified data).
12-25 Options are not presently used (March, 1973).
26-40 Reserved for individual plot options (not available).
GRID DIVISION AND GEOMETRY MODIFICATION
Application of the node concept based on input of grid division data
has been explained in the program description. Given (input) dimen-
sions of the aquifer are automatically adjusted by the program so
that the aquifer geometry is a multiple of the primary or coarse grid
size.
For BC1, the depth from the water table to the horizontal centerline
of the drain, H, is used for the coarse grid size. Both the depth
to barrier, DG, and drain spacing, AG, are modified. Since the user
has no control over the coarse grid size, the effect of geometry
adjustments is of interest. The computed drain discharge is a prime
indicator of this effect. Figure 24 presents the effect on the
dimensionless flow parameter of various dimensionless geometry
parameters:
a. Ratio of drain depth to barrier depth.
b. Ratio of drain spacing to drain depth.
c. Ratio of effective drain radius to drain depth.
For normal drainage situations, the solution is seen to be relatively
insensitive to changes in the aquifer geometry.
For BC2, the nominal size of the primary grid, HG, is input. Because
drain discharge is directly proportional to the drain spacing, A,
HG is adjusted so A is a multiple of the adjusted size. The barrier
depth, DG, is then modified. If HG is judiciously chosen, or if it
is made small relative to DG, adjustments can be minimized.
LIMITATIONS AND ERROR MESSAGES
The program description lists limitations on the size of the problem
that can be analyzed. They are repeated here:
169
-------
o
4.U
3.8
3.6
3.4
3.2
3.0
2.8
2.6
2.4
2.2
2.0
IB
1.6
1.4
1.2
i.O
^80, It
No
^40
V-20 :
Mo
x40-8
^-20
/4Q,e
^20
>0
,80,160
f===a.
3O,I60
- —
V|Q
^v.
^
_r
c
=====
r
d
=====
r
d
X
- - 0.40
W
"*-N
- = 0.20
=^— H
-=0.10
0,l60-[4- parameter
KIQ
T— -
_r
c
: = 0.05
\
\
"~^
- — ~^
V
X
^
'^-^,
water tabie-^
o '
Q=fl
.p
Denote
is 90
.deep
\
\
\
^
- . '• ' • -o
,,o • . •
Q
ow intc
er unit
s locat
% of th
aquifer
\
^/Cutc
t/ the
the
[£a_
) one drain
drain length
on (-§•) at whic
at for an infir
^ Kd
n'teiy
)ff limits based on
drain just touching
impermeable barrier
¥ °r |=%
*wwsnwMw*''W* A i
or f = (
\
\
x
^
\
ro
ro
- oo
\
\
en
o
V
\ [^
(Ti
D O.I 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
d/D
i \
l*r/d '
.0
FIGURE 24. DIMENSiONLESS FLOW
CIRCULAR TILE DRAINS BELOW HORIZONTAL WATER TABLE
170
-------
a. Maximum total number of nodes is 3,000 and the maximum number
in any vertical grid line is 400. Since these numbers depend on
the aquifer geometry and grid subdivision data, the user has
direct control. Equations to check these numbers prior to a
run are given in the next section.
b. Maximum number of interpolated points defining the streamlines
is 3,000. The user has indirect control on this number since it
depends on the number of nodes and streamlines, as well as the
aquifer geometry and boundary conditions. If this situation is
encountered, a coarser grid can be used or the number of stream-
lines reduced.
c. Maximum number of stream tubes is 20. If more are specified,
the program notes the condition, reduces the number to 20, and
continues. It should be noted that other programs in the water
quality model can only accept information for 10 tubes (i.e.,
SATCHEM and DEO200) .
d. Maximum seepage rate for BC2 must be less than the permeability.
This restriction is imposed by the physical nature of the problem
since the ultimate steady state would be the ponded water case.
Error messages are generated if the above limitations are violated.
In addition, failure of the infinite series for the solution constants,
stream functions, and potential functions produce error messages. A
list of all error messages is contained in the Appendix. Suggested
courses of action when errors do occur are included.
EQUATIONS FOR NUMBER OF NODES
If the user has specified subdivision information which causes more
than 3,000 nodes to be computed, an error message is generated and
execution is terminated. Should there be any suspicion that this
situation might prevail, the equations given below may be used to
avoid possibly unproductive runs.
For BC1, the equation is
NT = MP+(MP)(NP)+NP+(MH-N1)(L-1) [50]
+ (M1)(N1)(L2-1)+3L(2L+1)
where NT is the total number of nodes (excludes origin);
MP is the number of primary grid elements in the
horizontal direction;
NP is the number of primary grid elements in the
vertical direction; and
Ml, Nl, L are as described in the Data Deck section.
171
-------
For BC2, the equation is
NT = MP+(MP)(NP)+NP+(M1+N1)(L-1) [51]
+ (Ml)(Nl)(L2-l)+L.(3L+2)
where all terms are defined as before.
The determination of MP and NP requires use of the methods employed
by the program in modifying the aquifer geometry. For BC1, the
following procedure is followed:
a. Compute DG/H and AG/(2H).
b. Then NP is the integer part of DG/H if the decimal part is <0.5.
MP is the integer part of DG/H increased by 1 if the decimal part is
>0.5.
c. Similarly MP is the integer part of AG/(2H) is the decimal part
is <0.5.
MP is the integer part of AG/(2H) increased by 1 if the decimal
part is >0.5.
As an example, if DG = 28., AG = 150., and H = 10., then DG/H = 2.8,
AG/(2H) =7.5 from which NP = 3 and MP = 7.
For BC2, the following procedure is used:
a. Compute A/(2HG).
b. Then MP is the integer part of A/(2HG) if the decimal part is
<0.5.
MP is the integer part of A/(2HG) increased by 1 if the decimal
part is >0.5.
c. The modified primary grid size is H = A/(2MP).
d. Compute DG/H.
e. Then NP is the integer part of DG/H if the decimal part is
<0.5.
NP is the integer part of DG/H increased by 1 if the decimal
part is >0.5.
As an example, if DG = 28., A = 150., and HG = 10., then A/(2HG) = 7.5,
MP = 7, and H - 150/14 - 10.71. Computing DG/H = 2.6, NP = 3.
172
-------
An additional restraint on program use is the limitation that the
maximum number of nodes on a vertical grid line must not exceed 400.
The following equations may be used to compute this number:
BC1 NC = CNP+1) + Nl(L-l) + 2L f52l
BC2 NC = (NP+1) + Nl(L-l) + L [53]
LOGICAL ERROR MESSAGES
The following error messages may be generated during execution. The
boundary condition to which a particular error message applies is
indicated by the BC column. If the same message applies to both
boundary conditions, both numbers are shown.
No. BC_ Message
1 1 Computation of the stream function by subroutine SF0113
was attempted for a point outside the aquifer boundaries.
1 2 Computation of the stream function by subroutine SF0101
was attempted for a point outside the aquifer boundaries.
2 1 Computation of the stream function by subroutine SF0113
was attempted for the origin.
2 2 Computation of the stream function by subroutine SF0101
was attempted for the origin.
3 1 The series for the stream function did not converge
according to the tests applied in subroutine SF0113.
The maximum limit of 30 terms was reached.
3 2 The series for the stream function did not converge
according to the tests applied in subroutine SF0101.
The maximum limit of 999 terms was reached.
4 1,2 The number of nodes computed by subroutine SF0111 exceeds
the dimensioned storage of 3,000.
5 1 The series for the constant lower case Q did not con-
verge according to the tests applied in subroutine
SF0114. The maximum limit of 30 terms was reached.
5 2 The series for the constant C did not converge according
to the tests applied in subroutine SF0102. The maximum
limit of 999 terms was reached.
173
-------
No. BC Message
6 1 The series for the potential function did not converge
according to the tests applied in subroutine SF0115.
The maximum limit of 30 terms was reached.
6 2 The series for the potential function did not converge
according to the tests applied in subroutine SF0103.
The maximum limit of 999 terms was reached.
7 1,2 The number of points in a vertical column exceeds the
dimensioned storage of 400.
8 1,2 The number of interpolated points exceeds the dimen-
sioned storage of 3,000.
9 2 The number of points locating the streamlines exceeds
the dimensioned storage of 3,000 after adding points at
the water table.
10 1,2 The number of points defining a streamline exceeds the
dimensioned storage of 3,000.
11 1 The series for the potential function did not converge
according to the tests applied in subroutine SF0115.
The maximum limit of 30 terms was reached. The sub-
routine to compute the potential function was called
by SF0120.
11 2 The series for the potential function did not converge
according to the tests applied in subroutine SF0103.
The maximum limit of 999 terms was reached. The sub-
routine to compute the potential function was called
by SF0120.
12 1,2 The number of points defining a streamline exceeds the
dimensioned storage of 3,000 after subroutine SF0120.
13 2 The series for the potential function did not converge
according to the tests applied in subroutine SF0103.
The maximum limit of 999 terms was reached. Subroutine
SF0103 was called by SF0122.
14 2 The ratio of the seepage rate to the permeability is
equal to or greater than 1.00. There is no steady-
state solution under these conditions. The case of
uniform seepage must eventually yield the ponded
water case.
174
-------
No. BC Message
15 !»2 The number of required streamlines exceeds the dimen-
sioned storage of 20. This number is being reduced to
20 by the program.
Messages 1, 2, and 15 are not fatal. The remaining messages are
considered fatal and execution is terminated. If any of these error
messages are encountered, the following tabulation suggests means to
correct the deficiencies.
No. Course of Action
1>2 These are nonfatal and are primarily debugging messages.
Unless program changes are made, the user should never
encounter them.
3,5,6,11,13 These fatal errors all involve failure of the various
infinite series to converge. Details of the tests are
contained in the program description. Two alternatives
exist:
a. Increase the maximum number of terms in the
given subroutine. It is set at the start of the
subroutine by the arithmetic statement MMAX = 30
or MMAX = 999. Run times will increase.
b. Increase the test value for the ratio test
in the given subroutine. It is set at the start
of the subroutine by the arithmetic statement
RATI0 = .000001. Both run times and solution
accuracy will decrease.
The appropriate subroutines are included in the messages
and are summarized as:
BC Stream Solution Potential No. of
No. function constant function^ terms
30
999
Messages 3, 5, 6 are generated in computing values for
nodes or interpolated points. Subroutine SF0120 in
message 11 modifies points within the drain and adds
one on the perimeter. Subroutine SF0122 in message 13
computes water table locations for BC2.
1
2
SF0113
SF0101
SF0114
SF0102
SF0115
SF0103
175
-------
No. Course of Action
4,7 These fatal errors involve the aquifer geometry and grid
information. A redefinition of the grid using a coarser
net is required. The equations in the main part of this
manual can be used as guides. (See Equations for
Number of Nodes.)
8,9,10,12 These fatal errors involve both the grid data and the
number of streamlines. Either the grid should be made
coarser or the number of streamlines reduced. Mes-
sage 8 is generated during interpolation for all stream-
lines before adjustments at the drain or water table are
made. Messages 9, 10, and 12 apply to a single stream-
line. Message 8 will normally be encountered before any
of the others unless only 1 stream tube is used. Mes-
sage 9 is generated for BC2 when points are added to the
water table. Message 12 results when points are deleted
and added at the drain.
14 This fatal error involves use of the wrong boundary
condition based on the ultimate steady state. For
normal drainage applications, this situation indicates
severe waterlogging problems.
15 If more than 20 streamlines are desired, program
dimensions for V0L(I), XST(I), and TRAVEL(I) must be
increased. In addition, the error message generator
in subroutine SF0112 must be changed.
176
-------
SECTION XV
USER'S MANUAL
SATURATED CHEMISTRY PROGRAM
INTRODUCTION
The Saturated Chemistry Program was originally designed as part of
the irrigation return flow quality model. Dutt and Shaffer(4) have
presented the basic theory and concepts used. The program has been
expanded to allow its use in return flow studies not utilizing the
Unsaturated Flow and Chemistry Programs. Also, provision has been
included for salt movement across the barrier (lower boundary), and
user inputed values for ion exchange coefficients and CO partial
pressures. 2
The two-dimensional saturated system drained by parallel rows of
circular drains is divided into a number of stream tubes. Discharge
through each tube is identical although the soil mass is not. Indi-
vidual tubes may also be subdivided into a number of elements, each
having the same soil mass. As shown in Figure 25, these elements are
in series with respect to the path of flow.
The system is defined by specifying the number of tubes and for each
tube the total pore or water volume, average volumetric moisture con-
tent, stream tube width at the water table, and the number of sub-
divisions. Each tube can be initially treated as chemically homo-
geneous or each element can have a unique soils analysis. Once the
system is operating, an initially homogeneous system will soon become
heterogeneous due to different soil volumes in each tube, chemical
reactions, and piston displacement effects of the deep percolation
water. When a restart deck is punched at the end of a run, current
soils data for every element in every tube reflect the heterogeneous
condition.
Inputs to the saturated system include the amounts and quality of
deep percolation water. This may be taped output from the unsaturated
chemistry portion of the return flow model or may be input by cards.
When cards are used, inputs may represent output from the unsaturated
zone (either observed, predicted, or guessed data) or surface inputs.
Surface inputs, in turn, would logically be used when the unsaturated
zone is nonexistent (for example, the ponded water case), when the
zone is bypassed because changes within it are minor in comparison to
those in the saturated system, or when the program is being used to
simulate the unsaturated zone. In this case, input data could be the
estimated deep percolation and the corresponding quality of the applied
177
-------
Water table
Mid- point flow line
/Tube I
Tile drain
Segment or
element no.
In general, max.(I).
FIGURE 25. SEGMENTED FLOW TUBE IN SATURATED REGION.
EXAMPLE OF 10 SUBDIVISIONS. (MAX (I) = 10)
-------
water or water of surface quality can be entered along with the onfarm
irrigation efficiency. In the later case, the program removes the
consumptive use from the input water and concentrates it accordingly.
Output from the program includes printed results showing the quality
of water discharging from each tube, punched restart chemistry decks,
and taped quality output. The taped output is intended for input to
the drain effluent prediction portion (program DE0200) of the water
quality model.
UNIT NUMBERS
Card input is accomplished by READ statements without a unit number
and listed output by means of PRINT statements. The common input
(CIT) and output (C0T) devices are understood to be used. Punched
output (restart deck) is accomplished by PUNCH commands. In addition,
the program references additional units:
Unit Unit description
designation or purpose
15 The unit on which the file containing deep percolation
quantities and qualities is mounted. This file is
written by the unsaturated chemistry portion of the
water quality model. When these data are input by
cards, they are transferred to unit 15 and the program
functions as if a saved tape existed. (See Tape Input.)
16 The unit on which the file containing the quality of
water discharging from each tube of the saturated system
is mounted. This file is intended for input to the
drain effluent prediction portion of the quality model.
(See Tape Output.)
17 An internal scratch unit on which the last year's data
are written. This allows repeated reading of these data
as a single file.
18 Output data written to this file in format suitable for
a read in a subsequent run as unit 15 (IWR must = 1).
20 This file contains 7 coded output giving the total num-
ber of slugs generated for each stream tube during the
run. Information is used as partial input to the drain
effluent prediction submodel.
179
-------
DATA DECK STRUCTURE
There are six groups of input cards from which the data deck for a
given run is composed, as shown in Figure 26. The format for each
group, except group 6 (restart deck), is illustrated in Figure 26.
Variables are defined below. Note that the card numbers refer to
card types and not the actual number of a given type.
CARD GROUP 1 - Title Card
Column
No.
Var
type
Var
name
Card 1 of 1
1-80
TITLE
CARD GROUP 2 - Control Card
Card 1 of 2
1- 5
I ICARD
6-10
11-15
I INFIL15
I LSFIL15
16-20
I NREPT
21-25
I INFIL16
Variable
description
An 80-character alphanumeric title used
to identify the printed output.
Flag to indicate the input medium for
the deep percolation inputs.
If ICARD = 1, then read input from
cards (group 5) .
If ICARD $ 1, then read input from
magnetic tape (unit 15).
The initial file number to be read on
unit 15. (See Tape Input.)
The last file number to be read on
unit 15. When deep percolation data
are read from cards (ICARD = 1),
LSFIL15 is the number of years of data.
(See Tape Input.)
The number of times the last file on
unit 15 (LSFIL15) is to be repeated as
input. This option is available for
both tape and card input. (See Tape
Input.)
The number of the initial file to be
written on unit 16. (See Tape Output.)
180
-------
NANAL=0
r
ONE YEARS DATA, REPEAT
FOR SUCCESSIVE YEARS
y Restart deck
Soils analysis
GROUP 5 DETAILS
\
IRERUN = I
/lubedata
Salt from
barrier data
/Deep perc.
data
\
\
MRERUN^I
FIGURE 26. DATA DECK STRUCTURE
SATURATED CHEMISTRY SUBMODEL
181
-------
Column
No.
26-30
31-35
36-40
41-45
46-50
I
R
Variable
description
ISAME A flag to denote whether each tube is
assumed homogeneous with respect to
soil chemistry or if each segment in
the tubes has a different analysis.
If ISAME = 1, then each tube is homo-
geneous and there will
be one group 4 card
for each tube (MAXTU
cards).
If ISAME j* 1, then there will be a
group 4 card for each
element in a tube (MAX(I)
cards for tube I). (See
group 3 cards below.)
IRERUN Flag to denote a restart run.
If IRERUN = 1, then this is a restart
run. A restart deck
(group 6) is used and
initial soils analyses
are deleted (group 4).
If IRERUN t 1, then this is an initial
run and soils data are
inputed (group 4).
IPUNCH Flag to signal punching of a restart
deck.
If IPUNCH = 1, punch a deck at the end
of the run.
If IPUNCH •£ 1, then do not punch.
MAXTU Number of stream tubes used.
CHI Convergence criterion for the ion
exchange routine. It is the number
of moles/liter which is considered
as insignificant. Suggested values
range from 10~2 to 10~6 with a recom-
mended value of 10"1*. .
182
-------
Column
No.
51-55
56-60
61-65
66-70
71-75
76-80
Var Var
type name
I IDUMP
IINJE
SBYPAS
IREK
IDN
IPC02
Variable
description
If = 1, the chemical composition of the
water contained in each tube
element is printed at the end
of the run and represents the
chemical status at this point
in time.
If = 0, this output is suppressed.
If = 1, the salt injection option is
activated. (Group 3 must be
present.) See Figure 27 for
details.
If = 0, option is bypassed (group 3
must be absent).
If = 1, the chemical reactions (salt)
within the program are
bypassed.
If = 0, the chemical reactions are
included.
If = 1, cation exchange coefficients
for Ca-Na and Ca-Mg ion
exchange are read from the
soils analyses data cards.
If = 0, default values within the model
are used.
If = 1, the denitrification subroutine
is utilized.
If = 0, denitrification subroutine is
bypassed.
If = 1, user supplied pC02 values used1
(see group 6).
If = 0, program relate: pC02 to soil
moisture content.
Card 2 of 2
1- 5
IWR
If = 1, output written to file Tape 18
in binary format suitable for
a subsequent read as file
Tape 15.
If = 0, no write to file Tape 18.
183,
-------
Mass f lux,
* *
* *
« « X X * X
X "
x x x x „ x *
X XX X
* * x x" „* « x
* *
Water of aqu i fer
qual i ty (Iower
tube)
Transition zone(AX)
Barr i er
yater of Barrier
qual ity (constant
with t ime)
Mass Fluxj
Mass Flux;
D
D.
-££1
rf-X
Bas i c equat i on
Approximate equation
(steady state)
(for any t ime step)
Where: Ci is the concentration of constituent i.
D is a mean diffusion coefficient.
FIGURE 27. SALT INJECTION FROM THE BARRIER
-------
CARD GROUP 3 - Deep Percolation Data
This group is only used when ICARD = 1. It includes the three types
of cards described below. Cards must be placed in chronological order
with each year's data preceded by a type 1 card. This card is fol-
lowed by the appropriate combination of card types 2 and 3 as indicated
in Figure 26. •(See NANAL below.) Checks are not made on the chrono-
logical order and it is possible to include more than one application
on a single day by using separate type 2 cards. The total number of
years of data is LSFIL15. (See Tape Input.)
Column
No.
Card 1 of 3
1- 5
Card 2 of 3
1- 5
6-10
11-15
Var
type
Var
name
NIRR
IM0N,
IDAY,
IYEAR
16-25
R
APP
26-30
NANAL
Variable
description
Number of water applications or deep
percolation inputs for the current
year.
Deep percolation volumes.
The month, date, and year on which the
deep percolation occurred. For exam-
ple, July 17, 1973, would be IM0N = 7,
IDAY = 17, IYEAR = 1973. The month
and day are converted to the Julian
day, which along with the year, is
written on unit 15 for information
purposes.
Amount of deep percolation for this
application (inches3/inch2 or inches).
A + value is an input to the saturated
system, while a - value denotes flow
from the saturated to unsaturated
system.
Flag used to denote the need for a new
water quality analysis.
If NANAL = 0, the previously read
analysis is to be used
with this application.
If NANAL t 0, a new analysis is to be
read and a type 3 card
will immediately fol-
low this type 2 card.
185
-------
Column
No.
31-46
47-56
Card 3 of 3
1- 5
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
CARD GROUP 4
Card 1 of 2
1-10
Var
A
R
R
R
R
R
R
R
R
R
R
- Data
R
Var
name
S0URCE
E
ANH4
AN03
CA
AMG
ANA
HC03
C03
CL
S04
for Salt
DIFK
Variable
description
A 16-character alphanumeric name which
can be used to identify the source of
the deep percolation water. For exam
pie, PR0JECT WATER or PRECIPITATI0N
or SN0WMELT. It is printed on output
for information purposes only.
Onfarm irrigation efficiency (percent
basis) .
Quality analysis of deep percolation
water.
Ammonium (NH^*) concentration (meq/1).
Nitrate (N03~) Concentration (meq/1) .
Calcium (Ca++) concentration (meq/1) .
Magnesium (Mg*+) concentration (meq/1) .
Sodium (Na+) concentration (meq/1) .
Bicarbonate (HCOs") concentration
(meq/1).
Carbonate (C03=) concentration (meq/1) .
Chloride (Cl~) concentration (meq/1).
Sulfate (50^=) concentration (meq/1) .
Movement (Injection) From Barrier
Diffusion coefficient (cm2/day) values
11-15
DIFDIS
range from about 1 to 20.
Diffusion distance (cm).
186
-------
Column
No.
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
56-60
61-65
Var
type
I
I
R
R
R
R
R
R
R
R
Var
name
MIN
M0UT
SHCA
SHMG
SHNA
SHC03
SHHC03
SHCL
SHS04
SHN03
66-70
R TRTME
Card 2 of 2
1-10 R CLGT(L)
11-20
R CLGT(L+1)
Variable
description
First element of outside tube con-
tacting barrier (as measured from
surface of water table).
Last element of outside tube con-
tacting barrier.
Concentration of Ca** in barrier water
(mg/1)•
Concentration of Mg*4" in barrier water
(mg/1).
Concentration of Ha* in barrier water
(mg/1).
Concentration of C03= in barrier water
(mg/1).
Concentration of HC03~ in barrier water
(mg/1).
Concentration of Cl~ in barrier water
(mg/1).
Concentration of S0^= in barrier water
(mg/1).
Concentration of N03~ in barrier water
(mg/1).
Travel time to displace one slug in tube
contacting barrier (days).
Length of contact of element L with
barrier (cm).
Length of contact of element L with
barrier (cm).
187
-------
Column
No.
Var
type
Variable
description
CARD GROUP 5 - Stream Tube Data
Card 1 of MAXTU
1-10 R TUVEN(I)
11-20
21-30
Pore volume of tube (feet3/foot of
drain). Where 1=1 corresponds to
the tube closest the drain and
I = MAXTU to the tube farthest
from the drain.
R TUWEN(I) Width of tube at the water table (feet).
31-40
R THETA(I)
MAX(I)
Volumetric moisture content of the tube
at saturation. It is the ratio of the
volume of water to the volume of soil
(dimensionless).
Number of elements or segments that the
tube is to be divided into. (Program
limits this value to 25 or less.)
Note: One card is needed for each tube.
CARD GROUP 6 - Soils Analysis Data
This group is only included for an initial run (IRERUN ? 1) • There
will be one card per tube if each tube is homogeneous (ISAME = 1) or
one card per element in all tubes (MAX(I) card for tube I) if each
tube is nonhomogeneous (ISAME $ 1). Cards should be in order, start-
ing with tube 1 nearest the drain and ending with tube MAXTU farthest
from the drain. Where there is one card per element, cards for a
given tube should start with element 1 nearest the water table and
end at element MAX(I) nearest the drain. (See Figure 25.) Note that
element boundaries are not necessarily horizontal. This means that
care must be taken when assigning soil chemistry data to each element.
Most of the chemistry data comes from a soil extract analysis.
Card 1 of 2
1-5
R
CA(J)
(J denotes the element number.)
Calcium (Ca++) concentration (meq/1).*
Meq/1 denotes miHiequivalents per liter.
188
-------
Column
No.
6-10
11-15
16-20
21-25
26-30
31-35
36-40
41-45
46-50
51-55
Var Var
type name
R AMG (J)
56-60
61-65
66-70
71-75
76-80
R
R
R
R
R
R
R
R
R
R
R
ANA(J)
ANH4(J)
CL(J)
S04(J)
HC03(J)
C03(J)
AN03(J)
EC(J)
R CAL(J)
R XX5(J)
R XTRCT(J)
BD(J)
AKCS(J)
AKCM'(J)
Variable
description
Magnesium (Nig**) concentration
(meq/1).
Sodium (Na4"1") concentration (meq/1).
Ammonium (NHi/) concentration (meq/1).
Chloride (Cl~) concentration (meq/1).
Sulfate (SO^3) concentration (meq/1).
Bicarbonate (HC03~) concentration
(meq/1).
Carbonate (C03=) concentration (meq/1)
Nitrate (NOa") concentration (meq/1).
Cation exchange capacity (meq/100 gm
of soil) .
Indicator as to whether soil is cal-
careous (lime present).
If Cal(J) = 1.0, soil is calcareous.
If CAL(J) = 0.0, soil is noncalcareous.
Gypsum (CaS04 • 2H20) content (meq/100 gm
of soil) . Residual gypsum content at
moisture percentage of extract.
Ratio of water to soil at which the
soil extract was made, expressed
as a decimal (gm water/gm soil) .
Bulk density of the soil (gm of soil/ cm3
of soil) .
Ca-Na exchange coefficient (if the
default values are not wanted) .
Ca-Mg exchange coefficient (if the
default values are not wanted) .
189
-------
Column Var Var Variable
No. type name description
Card 2 of 2
1-10 R PC02(J) C02 partial pressure in atmospheres.
Restart Deck
These cards are punched at the end of a run if IPUNCH = 1. They are
reread on a restart run (IRERUN = 1) in which case initial soils data
(group 4) are deleted. The cards are punched by the statement:
PUNCH (IFMT) JYEAR, KC0UNT, (ANH4(J), AN03(J), UREA(J), CA(J),
ANA(J), AMG(J), HC03(J), CL(J), C03(J), S04(J), EC(J), XX5(J),
CAL(J), E5(J), C5(J), SA5(J), CAS0(J), AGS0(J), BNH4(J),
CMH201(J), BD(J), HK(J), AZE(J), J=l ,N) .
where IFMT = 2I5/(6E13.5/6E13.5/6E13.5/3E13.5, I5.E13.5)
followed by PUNCH(JFMT) (H0LD (L),L=1,10)
where JFMT = 6E13.5/4E13.5.
The above cards are repeated for each tube. It is seen that there are
3+4xMAX(I) cards per tube for a total of MAXTU (3+4xMAX(I)) cards for
all tubes. The variable definitions follow:
JYEAR - Current year number.
KC0UNT - Number of slugs out this tube so far.
J - Element or segment number within the tube. It goes from
1 to N where N is MAX(I).
ANH4(J) - Mass of ammonium nitrogen (NH^-N) in element (ug)*.
AN03(J) - Mass of nitrate nitrogen (N03-N) in element (ug).
UREA(J) - Mass of urea nitrogen (Urea-N) in element (ug).
CA(J) - Mass of calcium (Ca"1"*) in element (ug).
ANA(J) - Mass of sodium (Na+) in element (ug).
AMG(J) - Mass of magnesium (Mg**) in element (ug).
* ug denotes 10~6 gms or micrograms.
190
-------
HC03(J) - Mass of bicarbonate (HC03-) in element (ug).
CL(J) - Mass of chloride (Cl~) in element (ug).
C03(J) - Mass of carbonate (C03=) in element (ug).
S04(J) - Mass of sulfate (SOU=) in element (ug).
EC(J) - Exchange capacity of element (eq/gm).
XX5(J) - Gypsum (CaSCV2H20) concentration in element (moles/gm).
CAL(J) - Lime indicator for element as used on input for initial
soils analysis.
If CAL(J) = 1.0, lime is present.
If CAL(J) = 0.0, lime is not present.
E5(J) - Exchangeable calcium (Ca"1"*) in element (moles/gm).
C5(J) - Exchangeable magnesium (Mg++) in element (moles/gm).
SA5(J) - Exchangeable sodium (Ma*) in element (raoles/gm).
CAS0(J) - Undissociated gypsum (CaSO^) in element (moles/1).
AGS0(J) - Undissociated magnesium sulfate (MgSO^) in element
(moles/1).
BNH4(J) - Exchangeable ammonium (NH^*) in element (moles/gm).
CMH201(J) - Water content in element (cm3).
BD(J) - Soil bulk density of element (gm/cm3).
IIK(J) - Lime indicator for element as used internally.
If IIK(J) = 1, lime is present.
}•'•
AZE(J) - Value proportional to intercept of the equation relating
pC02 and soil moisture content.(5) A constant for each
soil computed from the extract analysis.
H0LD(L) - Accumulated masses of various inputs to the saturated
system that were left over after the last slug output
was computed. For future runs, these partial sums are
used to initialize the mass counters until a volume of
water sufficient to displace the water in one element
has been accumulated.
191
-------
H0LD(1) - Volume of water (cm3/cm of drain).
H0LD(2) - Mass of calcium (Ca"1"1") (ug/cm of drain).
H0LD(3) - Mass of magnesium (Mg++) (ug/cm of drain).
H0LD(4) - Mass of sodium (Na*) (ug/cm of drain).
H0LD(5) - Mass of ammonium nitrogen (NH^-N) (ug/cm of drain).
H0LD(6) - Mass of chloride (Cl~) (ug/cm of drain).
H0LD(7) - Mass of sulfate (S01+=) (ug/cm of drain).
H0LD(8) - Mass of bicarbonate (HC03~) (ug/cm of drain).
H0LD(9) - Mass of carbonate (C03=) (ug/cm of drain).
H0LD(10) Mass of nitrate-nitrogen (N03-N) (ug/cm of drain).
TAPE INPUT
Magnetic tape output from the unsaturated chemistry program (unit 2)
can be used to input deep percolation data to the saturated chemistry
program (unit 15). This tape contains the volume of water leached from
the unsaturated zone and the corresponding mass of salts leached by
constituents. Values are expressed on a per-unit area basis. The
tape is written without format (binary) and is structured as shown in
Figure 28.
Each file represents a year of data. In executing the saturated
chemistry program, only full years are used. To allow for restarts,
the user specifies the initial file (INFIL15) to be used as input for
the first year computed by this run as well as the last file (LSFIL15).
The program skips the first (INFIL15-1) files and positions for the
initial read of the first logical record in INFIL15. The program then
reads and executes through subsequent files until LSFIL15 has been
read.
This last file may be repeatedly reread for future years' input by
using the NREPT option. If NREPT - 1, the program rereads and uses
data in LSFIL15 a total of NREPT times. Logically, the user may wish
to use this option if the data in LSFIL15 represent either an average
year or the results when a "dynamic" equilibrium condition is attained.
When deep percolation data are input by cards (ICARD = 1), the inputs
are written on unit 15 using the same file structure as for magnetic
192
-------
B0T (beginning of tape
,Ioad po i nt
E0F (end of file
mar k
^
first
i
second
1 og i ca
last
records
^
»
'/.
i\
L.
^ '
^ !
1 f iIe/year, 1st year
2nd year
3rd year etc.
Each Iog i caI
yRITE(2)
where 1 YEAR
1 DAY
SUM0UT
DELN
AMT(I)
DEL(I)
record is written by the binary write statement:
1 YEAR, I DAY, SUM0UT, DELN, (AMT(I), DEL(I) I=I,IO)
is the year number
i s the Ju I i an day , ->
is the accumulative leachate volume (cnr/crrr)
is the deIta-leachate volume since previous
record was written (cm /cm )
is the cumulative mass of constituent 1 leached
since the start of the run per unit area (10 *
err/ or ug/gm )
is the mass of constituent I leached since the
last logical record was written per unit area
(10 °gm/cm or ug/crrr)
gm/
Const i tuents: I
1
2
3
4
5
Item I
nitrate-nitrogen (NOj-N) fe
ammonia-nitrogen (NH^-N) 7
urea-nitrogen (Urea-N) 8
calcium (Ca++) 9
sod i urn (No ) 10
I tern
magnesium
bicarbonate (HCO-p
chloride (Cl~)
carbonate (C0-i=)
sulfate (SO/, = )
FIGURE 28 ORGANIZATION OF TAPE 2 OUTPUT FROM UNSATURATED
CHEMISTRY SUBMODEL
193
-------
tape. In this case, inputs do not include urea. Individual card val-
ues correspond to the delta amounts written on the tape. Since these
are the only values used by the saturated chemistry program, the accu-
mulative values are not computed. The month and date are first con-
verted to the Julian day which is then written along with the year and
delta amounts on unit 15. Cumulative values plus the delta amount for
urea are written as dummy variables.
When card input is used, LSFIL15 is also the total number of files or
years of input data. The first year is written as file 1, the second
as file 2, etc. Once data are transferred to unit 15, INFIL15 and
NREPT are used as they were with tape input. It is also noted that
the user can equip a magnetic tape on unit 15 and save the created
files for future use.
Finally, it is noted that there are no restrictions on the number of
logical records in a file on unit 15. Although the unsaturated chem-
istry program outputs results on a daily basis, this is not required.
The program reads all records in a file, using an end-of-file mark
check to determine when a file has been exhausted. The number of
records in each file may also vary.
TAPE OUTPUT
Magnetic tape output is written on unit 16 and is intended for input
to the drain effluent portion of the water quality model. Individual
runs of SATCHEM produce one file on unit 16. This file contains
results starting with the first slug reaching the drain in tube 1 for
the present run, followed by the second slug in tube 1, continuing
until the last slug for tube 1 has been computed. Results for succes-
sive tubes up to MAXTU follow. When all tubes and slugs are finished,
an end of file mark (E0F) is written. Thus, individual files can
represent results for any number of years and the number of logical
records per file may vary between files.
Figure 29 illustrates organization of the files and logical records on
unit 16. Each logical record contains the tube number, slug number,
and the concentration of the constituents in mgm/1. The slug numbers
begin with slug 1 for the very first slug computed by the initial run
and run consecutively, even for restart runs.
To allow restarts, the user must specify the initial file to be writ-
ten on unit 16, INFIL16. The program then skips (INFIL16-1) files
and positions the tape for the first write on INFIL16. For the ini-
tial run, INFIL16 should be left blank or set to 1.
194
-------
301 (beginning of tape,
Ioad po i nt)
E0F (end of file mark)
\
'//
1 -\ \ 1 ? ?
1 — 1 1 1 i$
slug 1 slug 2
tube 1 tube 2
"It
^~"~-^ og
ast
slug
tube
slug
tube 2
-------
LIMITATIONS
Dimensioned storage in the program limits the problem size to: a max-
imum of 10 stream tubes; a maximum of 25 elements or subdivisions per
tube.
There are no internal program checks to assure that these limits are
not exceeded. In addition, only full years or complete files on
unit 15 are processed. It is noted that there is no limit on the
number of deep percolation inputs by cards or tape.
OUTPUT
Program output consists of printed listings, punched cards, and com-
puted results written on magnetic tape.
Printed data listed by the program on the common print file is stand-
ard output. It includes:
a. Deep percolation volumes and qualities if input is by cards
CICARD = i)
b. Control card data
c. Salt input from barrier data (IINJE = 1)
d. Constituent concentrations for each tube element at end of
run (IDUMP = 1)
e. Stream tube inputs
f. Initial soil analyses (if IRERUN ^ 1)
g. Outflow concentrations by constituents for each slug of every
tube
Punched card output is optional, being selected by IPUNCH =1. It
consists of the restart deck which may be reus'ed without modifica-
tion for reruns.
Taped output on unit 16 is standard. If a saved magnetic tape is not
desired, a physical tape should not be equipped on this unit.
196
-------
SIMULATION OF UNSATURATED ZONE
The Saturated Chemistry Submodel can be utilized to simulate both
chemistry and waterflow in the unsaturated zone. In cases where this
technique can be applied, the savings in data collection/preparation
and computer time can be substantial.
Coarse-textured and many medium-textured soils lend themselves to
simulation by this method. A comparison of results obtained using
the Unsaturated Flow and Unsaturated Chemistry Submodels with those
obtained using the Saturated Chemistry Submodel provides a means of
testing the efficacy of this technique for a particular soil.
In applying the submodel to the unsaturated zone, a single stream
tube can be used to represent a one-dimensional flow path from the
soil surface to the water table. The tube may be subdivided into
some fixed number of elements or segments with equal volumes and
spacing, or a single homogeneous segment can be selected.
Consumptive use of water within the root zone is assumed to occur in
the upper most segment and water movement is assumed to occur at some
characteristic "leaching" moisture content. The soil water content
at field capacity is often used for this purpose.
The use of multiple segments in the soil profile adds desirable numer-
ical dispersion to the simulation. About 10 segments per profile have
been found to add sufficient dispersion.
Output is written on tape 18 in a format suitable for a direct read
by the Saturated Chemistry Submodel as tape 15. In the latter case,
the submodel would be used to simulate the saturated zone.
The following discussion summarizes details necessary to run the Sat-
urated Chemistry Submodel in the unsaturated mode.
a. Input data can be read from tape 15 or from cards. If tape 15
is used, the consumptive use must already have been removed from the
applied water. The program can remove the consumptive use if cards
are utilized.
b. It is possible to assign a different initial chemistry to each
segment by setting ISAME t 1.
c. Set MAXTU = 1, unless more than one soil profile is being sim-
ulated per run.
197
-------
d. If IDUMP = 1, a "picture" of the unsaturated profile chemistry
will be printed at the end of the run.
e. Set IWR = 1 so that output file tape 18 will be written.
f. TUVEN(l) must be set equal to the occupied pore volume of the
flow tube taken vertically from the soil surface down to a fixed
water table.
For convenience, assume a 1-ft2 column of soil and compute the pore
volume occupied by water at field capacity or some other "leaching"
water content as follows:
TUVEN(l) = DEPTH * THETA(1V
fc
Where: DEPTH is the vertical depth to the water table (ft),
and THETA(1)_ is the volumetric soil water content at
fc
the "leaching" water content level.
g. Set TUWEN(l) equal to 1 ft.
h. Set THETA(l) equal to the value discussed under f. above.
i. Set MAX(l) equal to the number of segments desired in the
in the tube (e.g. - 10).
j. Data for the soil chemical and physical parameters should
be inserted into the input data deck with values for parameters
closest to the soil surface associated with the first element or
segment, etc. Values for PC02 and other parameters should be
inserted to reflect the unsaturated soil profile.
198
-------
SECTION XVI
USER'S MANUAL
DRAIN EFFLUENT PREDICTION PROGRAM
INTRODUCTION
This program is designed to combine or mix the water qualities and flow
discharges from individual stream tubes extending from the water table
to subsurface or surface drains.
Outputs from the drainout, saturated flow, and saturated chemistry pro-
grams are used as direct inputs to compute the quality and quantity of
the drain effluent. When the saturated chemistry portion is bypassed,
inputs from the unsaturated chemistry program are used.
In essence, the drain effluent program performs the routing process,
modifies steady-state results when monthly discharge data are used,
performs certain computations to obtain salt loads (masses), and gen-
erates the required output. This output includes listings and punched
cards.
PROGRAM APPLICATION
There are two computational methods included in the drain effluent
program. The method used depends on the manner in which the quality
input data are developed. Briefly, in one case the inputs represent
the quality of leachate water and in the other the quality of the dis-
charges from the individual flow tubes.
Output from either method includes printed results. An optional card
output is also available when monthly drain discharges are used. Con-
centrations of each parameter and salt load results can be selected.
There are many similarities in the inputs required by either computa-
tional method. However, understanding the basic differences between
the two types of quality inputs will be helpful in applying the drain
effluent program. The two cases are discussed below.
QUALITY OF LEACHATE WATER IS INPUT
When this approach is used, chemical transformations in the saturated
zone are not simulated and the saturated chemistry program is bypassed.
The drain effluent program merely routes the leachate water and its
salts down each stream tube to the drain.
199
-------
Steady-state flow is assumed and piston displacement theory is applied
to each tube. When the quality of the leachate water is changed at a
given time, a "front" dividing this water from that previously enter-
ing the saturated system is assumed to form. The front is assumed to
travel down each stream tube in accordance with the travel times com-
puted by the saturated flow program. The relative arrival times of
the front at the drain for each tube are then used to combine results
from each tube yielding the drain effluent quality.
Leachate water quality inputs can be entered by card input or by an
edited form of the taped output from the unsaturated chemistry pro-
gram. Although card input is usually limited to testing and "debug-
ing," they can be used to input observed field data for verification
purposes.
Because the final results are based on a steady-state (average) dis-
charge, outputs from the unsaturated chemistry program must be edited
so that each change in quality corresponds to the proper volume of
leachate water. Since the unsaturated chemistry program produces
both variable leachate volumes and qualities on a daily basis, it is
necessary to edit the data so that inputs are based on equal volumes
of water. The corresponding average concentrations of each parameter
are then computed. This is accomplished by a utility program called
EDITC0 (see Suggested Improvements).
When cards are used, unequal time steps may be used, with inputs only
being required each time the leachate water quality changes. However,
the volume of water is assumed equal to that which would exist under
steady-state drain discharge.
DATA FROM SATURATED CHEMISTRY PROGRAM IS INPUT
When the saturated chemistry program is used, each stream tube is sub-
divided into a number of equal volumed segments. Leachate water is
accumulated until the volume corresponds to the pore volume of the seg-
ment. The mass of each parameter is also accumulated, resulting in a
water having the average concentration during the given time period.
This volume of water with its associated quality is referred to as a
slug.
The saturated chemistry program then passes the slug down the tube
through successive segments, simulating the various chemical trans-
formations. The final result is the quality of water in each slug
which discharges into the drain from a given tube. It is this dis-
charging water quality that is placed on tape by the saturated chem-
istry program and used as input to the drain effluent program.
200
-------
Because chemistry changes in the saturated system are not simulated
when leachate water is input to the drain effluent program the dis-
at:rwa^rWat?hifality V*1"11* "«tlcal to tS.S'S'thJ lach-
ate water. This is true for every stream tube. However, when the
saturated chemistry inputs are used, each tube will have different
at least until a "dynamic"
The relative timing of each slug arrival is based on an average
drain discharge and the corresponding travel times. When this is
done, results represent average qualities for equal increments of
discharge. However, monthly drainage volumes determined by the
drainout program can be used to modify the steady-state results.
The variable discharge pattern is then used to compute the average
monthly concentrations.
UNIT NUMBERS
Six logical unit numbers are referenced when the quality of leachate
water is input: S, 6, 7, 8, 16 and 29. When the saturated chemistry
results are used as input there are 7+NST units referenced: S, 6, 8,
9, 16, 29, 62, and MUNIT = 16+1 to 16+NST where NST is the number of'
stream tubes. Since the maximum number of stream tubes is 10, the
maximum number of units is 17.
Unit numbers are:
Unit
designation
Unit description
or purpose
The common input tape or unit (CIT).
puter systems standard input device.
The common output tape or unit (C0T).
puter systems standard output device.
It is the corn-
It is the com-
The unit on which the "formula" or equation giving the
composition of the drain effluent is written (see
Suggested Improvements). It is only used when the
leachate water quality is input. This is normally a
scratch unit and is written with unformated (binary)
writes. However, it may be equipped and saved. Sub-
sequent runs must specify a saved tape by setting
I0C(7) = 1, 2, or 3.
The unit on which all data necessary for plotting con-
centrations or for computing the salt load results is
written. It is written in binary form and is usually
201
-------
Unit
designation
16
MUNIT
Unit description
or purpose
considered a scratch unit although it may be equipped
and saved for future use. It is used when input is
the leachate water quality or the drain effluent qual-
ity, although the structure and data are slightly dif-
ferent . For:
input of leachate water quality (I0C(4) = 1 or 2),
unit 8 is written when salt load results are
requested (I0C(5) = 1 or 2) or when a saved tape
is requested (I0C(6) = 1), When I0C(7) = 2 or 3,
then a previously generated saved tape on unit 8
is indicated and must be included in the run.
input of the drain effluent quality (I0C(4) = 3),
unit 8 is written when monthly concentrations are
to be computed (I0C(10) = 1) provided a request is
specifically made to create a saved tape (I0C(6) =
1). It is used for subsequent runs by setting
I0C(8) = 1.
The scratch unit on which the computed drain effluent
quality (composite of all tubes) for all parameters
is written each time a slug reaches the drain. It is
only used when the saturated chemistry results are
used as input (I0C(4) - 3).
The unit on which the tape containing the quality
inputs is mounted. Either leachate water is input
(I0C(4) = 2) or the saturated chemistry tape is
used (I0C(4) = 3).
Data on the magnetic tape mounted on unit 16 is not
formated and is in binary form. The tape structure
is shown in Figure 29 for the case where the quality
of leachate water is input and in Figure 30 when the
drain effluent quality is input.
The logical unit on which the water quality inputs
from the saturated chemistry program (as read from
unit 16 when I0C(4) = 3) are transferred. It is an
internal scratch unit on which the data is written in
binary form.
202
-------
NEXT PAGE
'PARAMETER UNITS
CODE
PARAMETER NAMES
(2nd. card if more than
6 parameters)
'PARAMETER NAMES
'NO. QUALITY
CHANGES B
PARAMETERS
Q ORIGINAL GW
UALITY (Znd.cardii
more than 8 param.
/ TRAVEL TIMES
/NO. STREAM
f TUBES
1
YES
-CARD TYPE
START OF INPUT DECK
FIGURE 30. DATA DECK STRUCTURE -DRAIN EFFLUENT PREDICTION SUBMODEL
203
-------
END OF DATA DECK
MONTHLY DISCHARGE
DATA {2 cords per
FIGURE so. (CONTINUED;
204
-------
Unit
designation
Unit description
or purpose
The data for each stream tube are placed on different
units with MUNIT starting one higher than 16, being
incremented by one for each additional tube. Thus:
MUNIT =16+1 for stream tube 1
MUNIT =16+2 for stream tube 2
MUNIT a 16 + NST for stream tube NST
Since NST has a maximum value of 10, MUNIT will range
from 17 through 26.
29 The unit to which plotter commands are written when
the plot subroutines are included in the model. This
corresponds to the default unit number when the CRT
unit is not specified by the user.
62 The card punch unit of the computer system.
DATA DECK STRUCTURE
There are 23 types of input cards in the input deck, as shown in Fig-
ure 30. The variables are defined below. It should be noted that not
all card groups are included in a given run.
CARD GROUP 1 - Title Information
Column
No.
Var
type
Card 1 of 1
1-80 A
TITLE
Variable
description
An 80-character alphanumeric title used
to identify all printed and plotted
output. If this title is desired to
be centered on output, it must be cen-
tered on the card.
205
-------
Column
No.
Var
type
Variable
description
CARD GROUP 2 - Problem Selection
Card 1 of 1
1 I
2 I
40
IPR0B(1) A control array which calls the proper
IPR0B(2) overlay for execution. Individual
programs in the water quality model
were originally designed for inclu-
sion in an overlay system. Presently
IPR0B(40) only the saturated flow and drain
effluent portions are included. When
IPR0B(I) is set to 1, the Ith overlay
is called. The program will accept
up to 40 entrys. Thus:
IPR0B(1) = 1, call the saturated flow
routine
IOR0B(2) = 1, call the drain effluent
routine
On output, values of IPR0B(I) are listed
as I0C values under the heading PROB-
LEM CONTROL INFORMATION FOR MAIN
PROGRAM.
CARD GROUP 5 - Input/Output Options
Card 1 of 1
1 I
2 I
I0CC2)
I0C(40)
An input/output option array used to
specify the type of inputs, the
required computations and the
printed and plotted output. (See
Input/Output Options for details.)
40 I
CARD GROUP 4 - Number of Stream Tubes
Card 1 of 1
1-10 I
NST
The number of stream tubes to be used.
When the saturated chemistry output
is used, NST must correspond to the
number used in the saturated chemis-
try program.
2O6
-------
Column
No.
Var
type
Var
name
Variable
description
CARD GROUP 5 - Travel Times
Card 1 of 1
1-8 R
9-16 R
73-80 R
TRAVEL(1)
TRAVEL(2)
TRAVEL(10)
Time of travel assuming piston displace-
ment down tube I (years). 1=1
corresponds to the tube nearest the
drain, I = 2 to the next closest tube,
etc. These values are computed by the
saturated flow program and are part of
its output. Up to 10 travel times can
be included.
CARD GROUP 6 - Number of Quality Changesand Parameters
Card 1 of 1
1-5
6-10
R
NUQUCH The number of changes in the quality of
the leachate water when either cards
are used to input quality (I0C(4) = 1)
or when input is from the unsaturated
chemistry program (I0C(4) = 2).
When input is the saturated chemistry
program results (I0C(4) « 3), NUQUCH
is the maximum number of slugs leav-
ing any tube. It will always corre-
spond to tube 1. The value is obtained
from the saturated chemistry programs
printed output.
NUQUPA The number of quality parameters being
considered. It includes the number
of anions plus the number of cations,
but does not include the total load.
NUQUPA must correspond to the number
of parameters used in the unsaturated
or saturated chemistry programs when
their output is used. NUQUPA must be
in the range 1 <_ NUQUPA £ 15.
207
-------
Column
No.
Var
type
CARD GROUP 7 - Parameter Names
Card 1 of N
1-10 A APANAM(l)
71-80 A APANAMC8)
(All fields may not be used)
Card N of N
1-10 A
Variable
description
A 10-character alphanumeric name used
to identify the quality parameter,
where I « 1 to NUQUPA. (N = 1.) If
NUQUPA exceeds 8, two cards are
required (N * 2): the first giving
APANAM(l) - APANAM(8) and the-second
giving APANAM(9) - APANAM(NUQUPA).
These labels are used to identify all
printed and plotted output, and should
be in the same order as the quality
inputs coming from cards (I0C(4) = 1)
or from the chemistry programs
(I0C(4) = 2 or 3).
APANAM(9)
61-70 A APANAM(IS)
(All fields may not be used)
CARD GROUP 8 - Parameter Units Code
Card 1 of 1
1 I
2 I
lUNIT(l)
IUNIT(2)
15 I IUNIT(15)
(All fields may not be used)
A numeric units code for parameter I
where the order is identical to that
for APANAM(I) and the quality inputs,
If IUNIT(I) -
IUNIT(I)
1, the units are in
milligrams/liter.
(On output they
are labeled MG/L)
the units are in
milliequivalents/liter.
(On output they are
labeled MEQ/L.)
208
-------
Column
No.
Var
type
Variable
description
CARD GROUP 9 - Original Ground Water Quality
Card 1 of N
11-18
19-26
R
R
QUAL(l)
QUAL(2)
67-74 R QUAL(8)
(All fields may not be used)
Card N of N
11-18 R
19-26 R
QUAL(9)
QUAL(IO)
59-66 R QUAL(15)
(All fields may not be used)
CARD GROUP 10 - Reference Time
Card 1 of 1
1-10 R
TIMREF
The initial or original quality of the
ground water at the start of the run
for parameter I where the order is
identical to that for APANAM(I) and
the other quality inputs. Units for
the initial quality must correspond
to those specified by IUNIT(I).
If NUAUPA <_ 8, N = 1
If NUAUPA > 8, N => 2
If there are more than 8 parameters,
2 cards must be used with the first
giving QUALI(l) - QUALI(8) and the
second QUALI(9) - QUALI(NUQUPA).
It is noted that the original ground
water quality is only used when the
saturated chemistry results are not
used. Since chemical transformations
are not simulated and piston displace-
ment is assumed, the leachate water
is merely routed to the drains. The
aquifer is assumed to be homogeneous
with respect to water quality and the
quality of the drain discharge will
be identical to the initial quality
until the first front reaches the
drain.
The reference time to be used as the
initial time at the start of the run
(years). It must be consistent with
times on the type 11 and type 12
cards. TIMREF could be 0.0, 1974,
etc.
209
-------
Column
No.
Var
type
Variable
description
CARD GROUP 11
Leachate water quality changes Group N cards must be arranged in chron-
ological order, starting with the first change and ending with the last
change (NUQUCH).
Card 1 of N
1-10
TIQUCH(I)
11-18
19-26
R
R
PPM (1,1)
PPM(I,2)
67-74 R PPM(I,8)
(All fields may not be used)
Card N of N
The time of water quality change I
(years). It must use the same time
base as TIMREF. Thus, if TIMREF =
1974 and the first change is half-way
through the first year, TIQUCH(l) =
1974.5 years.
The quality of parameter J for the Ith
change. Units correspond to those
specified on the group 8 card. Param-
eters are arranged in the same order
as APANAM(J) and the other inputs.
It is noted that PPM(I,J) is the
actual quality at this time, rather
than an incremental or delta value.
If there is no change for a given param-
eter at this time, the corresponding
PPM(I,J) may be left blank. The pro-
gram will then use the previous value.
Thus, only parameters actually experi-
encing changes need be punched on
cards.
If there are more than 8 parameters,
two Group 11 cards are required. The
first card will contain quality for
parameters 1 through 8 while the sec-
ond contains quality for parameters
9 through NUQUPA.
210
-------
Var
Column
No.
CARD GROUP 12 - Finishing Time
Card 1 of 1
Variable
description
1-10
R
TIMFIN
The time at which the run is terminated
(years). It must use the same time
base as TIMREF. For example, if a
16-year run is desired and TIMREF =
1974, TIMFIN = 1974 + 16. = 1990.
CARD GROUP 13 - Unsaturated Chemistry Tape Information
Card 1 of 1
1-5
NINTYR The number of equal time intervals
per year on the special unsatur-
ated chemistry tape produced by
the utility program EDIT0. (See
Figure 30 for tape for mat.)
CARD GROUP 14 - Saturated Chemistry Tape Information
Card 1 of 1
1- 5 I
NINTYR
6-10
NFILE
The number of elements or segments that
each stream tube was divided into in
running the saturated chemistry pro-
gram. (For example, if 0.1 pore vol-
umes were used, NINTYR would be 10.)
It is noted that although the satu-
rated chemistry program allows for a
different number for each tube, the
drain effluent program requires the
same number.
The number of files of data on the sat-
urated chemistry output tape that are
to be read [3, Tape Output]. Individ-
ual runs of program SATCHEM produce
one file. The organization of each
file is shown in Figure 29.
211
-------
Column
No.
Var
type
Variable
description
CARD GROUP 15 - Number of Slugs for Each Tube on Tape
Card 1 of 1
1- 5
6-10
N0SLUG(1)
N0SLUG(2)
46-50 I N0SLUG(10)
(All fields may not be used)
The total number of slugs of water
reaching the drain that are contained
on the saturated chemistry input tape
where 1=1 corresponds to the tube
adjacent to the drain, I = 2 to the
next closest, etc. The values of
N0SLUG(I) are obtained from the
printed output of the saturated chem-
istry program, or from tape unit
written by saturated chemistry.
CARD GROUP 16 - Combining Weights
Card 1 of N
1-10
11-20
R
R
C0MWT(1)
C0MWT(2)
71-80 R C0MWT(8)
(All fields may not be used)
Card N of N
1-10 R
11-20 R
C0MWT(9)
C0MWT(10)
61-70 R C0MWT(15)
(All fields may not be used)
The combining weight of parameter I
where the order is identical to that
for APANAM(I) and the other quality
inputs. The combining weights are
equal to the atomic weight divided
by the valence and are used to con-
vert qualities in milliequivalents/
liter to milligrams/liter.
When quality inputs are by cards and
the units specified by IUNIT(I) are
in milliequalvalents/liter for one
or more parameters, then the combin-
ing weight card is required. Only
the combining weights for those
parameters expressed in MEQ/L need
be included with the others left
blank.
When there are more than 8 parameters,
two cards must be included. The
first will contain combining weights
for parameters 1-8 and the second
for parameters 9 - NUQUPA.
212
-------
Column
No.
Var
type
Variable
description
CARD GROUP 17 - Steady-state Drainage Data
Card 1 of 1
1-10
11-20
21-30
R Q The steady-state drain discharge
(cubic feet/year/foot of drain). It
is the total discharge from the aqui-
fer on both sides of the drain. The
numerical value can be obtained from
the saturated flow program's printed
output.
R SPACE The horizontal drain spacing (feet).
R V0LUME The volume of water stored in the aqui-
fer from the water table to the bar-
rier between drains (cubic feet/foot
of drain). When all salt load com-
putations are requested (I0C(5) = 2),
the V0LUME must be included. Other-
wise the program will not use it.
CARD GROUP 18 - Total Load Array
Card 1 of 1
1 I
2 I
ISARRY(l)
ISARRY(2)
15 I ISARRY(15)
(All fields may not be used)
A flag used to indicate which parameters
are to be included in the T0TAL L0AD,
where I is the parameter in the same
order as APANAM(I) and the other qual-
ity inputs. If:
ISARRY(I) = 1, then include this param-
eter in the T0TAL L0AD.
ISARRY (I) ? 1, then exclude this param-
eter from the T0TAL
L0AD.
CARD GROUP 19 - Number of Years of Monthly Discharge Data
Card 1 of 1
1- 5 I
NYEARS The number of years of monthly discharge
data that are to be read. These data
are contained on the type 21 cards
which are output from the drainout
program.
213
-------
Column Var Var Variable
No. type name description
CARD GROUP 20 - Steady-state Discharge
Card 1 of 1
1-10 R Q The average yearly drain discharge used
to calculate arrival times in the sat-
urated flow program (acre-feet/acre/
year).
CARD GROUP 21 - Monthly Discharge Data
These cards are output from the drainout program, except for the title
card which is eliminated.
Card 1 of M
1-5 I IYEAR The year number of the data.
6-7 I ICRD The card order within the year. ICRD = 1
contains data for the first 6 months
and ICRD = 2 contains data for the
second 6 months.
8-19 R QM0N(1) The monthly drain discharge for month I
20-31 R QM0N(2) where I = 1 is January, I = 2 is Feb-
ruary, etc. (acre-feet/acre/month).
• * *
* * *
R
CARD GROUP 22 - Microfilm Identification
The model version included in Volume V does not include plot sub-
routines. However, a blank card must be added here.
CARD GROUP 23 - Fiche Header Card Information
Fiche capability not included in this version of the model. No card(s)
needed.
Input/Output Options
Individual input, computation, and output options are specified by
setting the appropriate I0C(I) variable to the desired value on the
214
-------
group 3 card. For many options, the variables act as switches with
a value of 1 turning the option on and any other value leaving it off.
Except for I0C(21) which controls card output individual printed output
options are selected by the I0C(11) - I0C(25) variables. An individual
option is selected by setting the I0C(I) value to 1. All of these
options, except for I0C(21) can be selected by setting 100(1) = 1.
This eliminates the need of turning the individual options on.
Plotted* output options are selected by the I0C(26) - I0C(40) variables.
The two basic methods of inputing water quality changes are selected by
the I0C(4) option. When I0C(4) = 1 or 2 the quality of leachate water
is used (i.e. water leaving the unsaturated zone) and if I0C(4) = 3 the
drain effluent quality is used (i.e., results computed by the saturated
chemistry program). Some options are used for both cases while others
apply to only one. Tables 16 and 17 summarize this information and
include a list of options that are presently unused. Restrictions in
the use of each option are also noted in the following list.
Option
Select all print options, I0C(11) - I0C(25).
(I0C£21) is not included since it is a punch
option)
2 1 Select all plot options, I0C(26) - I0C(40).
3 1 Read all required saturated flow data from cards.
This includes the number of stream tubes (card
type 4), the travel times (card type 5) as well
as the steady-state discharge, drain spacing and
volume of water in storage (card type 17) when
the leachate water quality is input (I0C(4) =
1 or 2) and salt load computations are requested
(I0C(5) = 1 or 2).
3 ^i The above data are expected to be in C0MM0N stor-
age as a result of running the saturated flow pro-
gram first.
4 _ This option controls the manner in which quality
inputs are made to the program.
4 i The quality of the leachate water is input via
cards (card type 11).
* Plot routines are available, but were not included in the model,
215
-------
TABLE 16. DRAIN EFFLUENT PREDICTION PROGRAM
Input/Output Option Summary
Available options
IOC(I)
option
1 =
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Leach ate
input
IOC (4) = 1 or 2
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
Drain effluent
input
IOC (4) = 3
X
X
X
X
X
X
X
X
X
X
X
X
X
X
Unused
options
X
X
X
X
Graphical
output options
not included in
this version of
the model.
216
-------
TABLE 17. DRAIN EFFLUENT PREDICTION PROGRAM
Salt Load Considerations When
Drain Effluent Quality is Input
(I0C(4) = 3)
Accumulative
Item
Compute
Print
Plot
Monthly
I0C(5) = 1
or
I0C(5) = 2
I0C(17) = 1
**
Annual
I0C(5) = 1
or
I0C(5) = 2
I0C(18) = 1
**
Rates*
Monthly
I0C(5) = 2
I0C(17) = 1*
**
Annual
I0C(5) = 2
I0C(18) = 1*
**
* Rates cannot be output unless computed. Thus, if
I0C(5) ? 2, the program ignores the I0C(17), I0C(18), and
I0C(32) commands to output nonexistent results.
** Plot capability not included in this version of the model,
but can be added (see Graphical Output).
Option
The quality of the leachate water is input via
magnetic tape from the special ability program
EDITC0 - mounted on unit 16, with the format
given on Figure 30.
The quality of the drain effluent is input via
magnetic tape from the saturated chemistry pro-
gram SATCHEM. This tape is mounted on unit 16
with the format as indicated in Figure 29.
Quality inputs are not inputed and results will
be meaningless.
In general, this option is used to select the
computation of salt load information. The spe-
cifics depend on whether quality inputs are for
the leachate water or for the drain effluent.
The I0C(5) options are now discussed separately
for each method of input.
any other
value
217
-------
Options
Leachate water quality is input (I0C(4) = 1 or 2).
Compute the accumulated salt load (tons/acre) in
the drain effluent plus the average discharge rate
of salt (tons/acre/year) on an annual basis. The
average annual rates for the entire run are also
computed. This is done for the individual param-
eters and the T0TAL L0AD.
Compute the accumulative salt load (tons/acre)
and the average discharge rates of salt (tons/
acre/year) for both the leachate water and the
drain effluent. The mean annual rate for the
entire run is also included. Results are for the
individual parameters and the T0TAL L0AD.
In addition, if printed results are selected
(I0C(17) =1), mass balance results are included
for the aquifer. The initial load in storage
(tons/acre) plus the change in storage and the
salt remaining in storage (tons/acre) are com-
puted for each year of the study on an accumu-
lative basis. This is done for each parameter
and the T0TAL L0AD.
any other This option is bypassed.
value
Drain effluent quality is input (I0C(4) = 3).
If a saved tape 8 from a previous run is not
being used (I0C(8) $ 1), the I0C(10) = 1 option
must be included for the run. Logically, both
the I0C(8) = 1 and I0C(10) = 1 options should
not be used in the same run. If they are, data
from tape 8 supercedes that computed by the
I0C(10) = 1 option for this run.
Compute the accumulative salt load in the drain
effluent (tons/acre) using monthly data on both
a monthly and on an annual basis. Results for
the individual parameters as well as the T0TAL
L0AD are included. Monthly results are listed
if I0C(17) = 1 while annual results are listed
218
-------
Options
if I0C(18) = 1. Both monthly and annual results
are plotted if I0C(30) =1.*
Compute the same results as for I0C(5) = 1 and
in addition compute the rate of salt load removal
in the drain effluent (tons/acre/year) using
monthly data on both a monthly and on an annual
basis. Results for the individual parameters,
as well as the T0TAL L0AD, are included. Monthly
results are listed if I0C(17) = 1 while annual
results are listed if I0C(18) = 1. Both monthly
and annual results are plotted if I0C(32) =1.*
any other This option is bypassed.
value
Write a saved tape 8 containing the date of the
run, title card, parameter name, parameter units
code, parameter units label, and the computed
concentration for each time and parameter. The
data are written in unformated form, with a sep-
arate file structure depending on the manner of
quality inputs. When the drain effluent quality
is input (I0C(4) = 3), the I0C(6) = 1 option is
only available when monthly drain discharge data
are used (I0C(10) = 1). There are no restrictions
when the leachate water quality is input (I0C(4) =
1 or 2).
The saved tape 8 may be read in subsequent runs by
proper use of the I0C(7) and I0C(8) options when
the quality of leachate water is input. When the
drain effluent quality is input, tape 8 usage is
controlled by the I0C(8) option.
This option is only available when the leachate
water quality is input (I0C(4) = 1 or 2). When
the drain effluent quality is input (I0C(4) = 3),
this option is bypassed regardless of the value
of I0C(7).
* When plots are requested, the accumulative results and rates appear
on the same graph. Because of the scaling and plotting process, the
salt load discharge rates cannot be plotted unless the accumulative
salt load is. Consequently, when I0C(5) = 2 and I0C(32) = 1, the
accumulative salt load is also plotted even though I0C(30) / 1.
219
-------
Value Options
An "equation" or "formula" for the quality of the
drain effluent in terms of each change in quality
of the leachate water is first developed.*
These compositions or results are written on
unit 7. Normally unit 7 is a scratch unit.
However, when there are many changes, data may
overflow the disks requiring that a tape be used.
This slows program execution because the unit is
rewound for each parameter.
1 Use a previously prepared tape 7. However, all
data must be read as though this tape were not
available, including quality changes. Computa-
tions start with the first parameter, the concen-
tration of each being computed with time. The
I0C(8) value should be set to zero or left blank
in this case.
2 Follow the procedure indicated by the I0C(7) = 1
option, except computations do not resume with
the first parameter. Instead, computations begin
with the parameter following the last one com-
puted and already written on a saved tape 8. The
number of the last parameter is given by the
value of I0C(8).
3 Skip the generation of information on tapes 7 and
8 and go directly to the computation of the salt
load results selected by the I0C(5) option. A
saved tape 8 (but not tape 7) is required. All
inputs, including quality changes, are still
required.
any other Follow the normal computation process starting
value at the beginning of the problem. Data will be
generated and written on unit 7.
? The use of the I0C(8) variable depends on the
manner of quality inputs. Where the leachate
water quality is input, the I0C(8) value should
be set to the number of quality parameters for
which results have been written on a saved tape 8,
The program reads over the required number of
* The present program method is computationally inefficient when the
number of leachate water quality changes is large. (See suggested
improvements.)
220
-------
Value
any other
value
records on tape 8 in preparation to adding new
data.* (See I0C(7) option.)
If I0C(7) * 0 or 1, I0C(8) should be zero or
blank If 100(7) . 2, then I0C(8) must be set
properly. If I0C(7) = any other value, the
value of I0C(8) has no effect.
When the drain effluent quality is input
UPH.4) = 3), a saved tape 8 is to be read.
Computations commence with salt load calcula-
tions. Thus,, the I0C(8) option is logically
used only when I0C(5) = 1 or 2. The input
deck is then truncated. Tape 8 mush have been
created in a previous run by setting I0C(10) = 1
and I0C(6) = 1 when I0C(4) = 3.
This option is not used.
any other
value
10
This option is only available when the leachate
water quality is input (I0C(4) = 1 or 2) and
salt load computations are requested (I0C(5) = 1
or 2).
Compute the average yearly concentrations of the
drain effluent when I0C(5) = 1 or 2 and/or the
leachate water when I0C(5) = 2. Results can
then be printed using I0C(19) and plotted using
the I0C(33) and I0C(34) options.
This option is not used.
This option is only available when the drain
effluent quality is input (I0C(4) = 3).
* Note that the 11 format allows only a single-digit number. There-
fore, if more than 9 parameters are on tape 8, the program must be
modified. (See Limitations and Suggested Improvements.)
221
-------
L Value Options
Read monthly drain discharge data (card type 21)
(i.e., output from drainout program) and com-
pute monthly concentrations. This option is only
used when the quality of the drain effluent is
input. Additional inputs include card types 18,
19, 20, and 21. Results are computed for each
parameter and the T0TAL L0AD. They are listed
when I0C(20) = 1, plotted when I0C(26) = 1
and/or I0C(27) = 1, and punched when I0C(21) = 1.
Results are placed on a saved tape 8 when
I0C(6) = 1.
10 any other This option is not used.
value
11 1 Print travel times as read from card type 5.
12 1 Print the leachate water quality for all changes
inputed (I0C(4) = 1 or 2). If I0C(4) = 3, this
option is ignored.
13 1 Print the arrival time of each piston displace-
ment front as it reaches the drain for each tube.
This is primarily a "debug" print. Only the
first 16 fronts are printed.
When the drain effluent quality is input
(I0C(4) = 3), the word front on output corre-
sponds to slug or the boundary between succes-
sive slugs. If there are less than 16 slugs/
fronts in the first tube, only NUQUCH are out-
puted. If there are less than 16 slugs/fronts
(or NUQUCH if it is less than 16) for any of the
other tubes, printed results may not be meaningful.
14 1 Print the composition of the drain effluent in
terms of the various leachate water quality
inputs (I0C(4) = 1 or 2) at the time of arrival
of each front. This is the information written
on unit 7. This is primarily a "debug" print.
Only the first 16 terms in the equation are
written. This option is not available when
I0C(4) = 3.
222
-------
Options
Print the values in the sorting arrays giving
the arrival times, corresponding front, and the
corresponding stream tube numbers. This is a
"debug" print and has no special headings. This
option is available when I0C(4) = 1 or 2, but
not for I0C(4) = 3.
16 1 Print the concentration of the drain effluent
for each parameter whenever a displacement front
reaches the drain. This option is only available
when I0C(4) » 1 or 2, but not for I0C(4) = 3.
!7 1 Print the load results that are computed by the
I0C(5) option when I0C(4) = 1 or 2. When
I0C(4) = 3, print the monthly load results that
are computed by the I0C(5) option.
18 1 This option is only available when I0C(4) = 3.
Print the annual load results that are computed
by the I0C(5) option.
19 1 Print the average yearly concentrations selected
by the I0C(9) option (i.e., drain effluent if
I0C(5) = 1 and drain effluent plus leachate water
if I0C(5) = 2). This option is only used when
I0C(4) = 1 or 2, but not I0C(4) = 3.
20 1 Print the monthly concentrations for each param-
eter plus the T0TAL L0AD when I0C(4) = 3 and
I0C(10) = 1. This option is not used when
I0C(4) - 1 or 2.
21 1 Punch the monthly concentrations for each param-
eter and the total load as well as the volume of
drain discharge. (See Card Output.) This option
is only available when I0C(4) = 3 and I0C(10) = 1.
22-25 - Not presently used.
26-34 _ Graphical output options; necessary subroutines
have not been included in this model (see Graphi-
cal Output).
35-40 - Not presently used.
223
-------
Graphical Output
A graphical output package is available for use with the model but has
not been included as part of Volume V. The package can be obtained
from the U.S. Bureau of Reclamation, Engineering and Research Center
in Denver, Colorado, Code 750.
224
-------
REFERENCES
(1) Dumm, L. D-, New Formula for Determining Depth and Spacing of
Subsurface Drains in Irrigated Lands. Presented at Annual
Meeting of ASAE, Chicago, Illinois, December 7-9, 1953.
(2) Dumm, L. D., Validity and Use of the Transient-Flow Concept
in Subsurface Drainage. Presented at 1960 Winter Meeting
of ASAS, Memphis, Tennessee, December 4-7, 1960. 28 p.
(3) Dumm, L. D., The Transient-Flow Theory and Its Use in Sub-
surface Drainage of Irrigated Land. Water Resources Confer-
ence, Irrigation and Drainage Division, ASCE, New York,
October 16-20, 1967.
(4) Dutt, G. R. and Shaffer, M. J., A Computer Model for Predicting
Nitrate and Other Solutes of Agricultural Drain Water. Final
Report to U.S. Bureau of Reclamation. The University of
Arizona, 1972.
(5) Dutt, G. R., Shaffer, M. J., and Moore, W. J., Computer
Simulation Model of Dynamic Bio-Physicochemical Processes in
Soils. Technical Bulletin 196. Agricultural Experiment
Station. The University of Arizona, 1972. 128 p.
(6) Glover, R. E., Ground-Water Movement. Engineering Monograph
No. 31. USDI. Bureau of Reclamation, 1964. 67 p.
(7) Jensen, M. E., Empirical Methods of Estimating or Predicting
Evapotranspiration Using Radiation. Evapotranspiration and
Its Role in Water Resources Management. ASCE Conference
Proceedings. December 5-6, 1966, p 49-53 and 64.
(8) Jensen, M. E., Wright, J. L., and Pratt, B. J., Estimating
Soil Moisture Depletion from Climate, Crop, and Soil Data.
Transactions, American Society of Agricultural Engineers.
September-October, 1971. p 954-959.
(9) Kunze, R. J., Uehara, G., and Graham, K., Factors Important
in the Calculation of Hydraulic Conductivity. SSSAP.
32:760-765, 1968.
(10) Land Drainage Techniques and Standards. Tentative. Office of
Chief Engineer, Bureau of Reclamation, Denver, Colorado,
February 1964.
225
-------
(11) Moody, W. T., Memorandum to Chief, Division of Drainage and
Ground-water Engineering, Subject: Table for Computation of
Hooghoudt's Equivalent Depth. July 17, 1964. 5 p.
(12) Ribbens, R. W., Report on Water Quality Model Developed in Con-
nection with Interagency Wastewater Studies, San Luis Drains,
California. USDI, Bureau of Reclamation, Engineering and
Research Center, Denver, Colorado, September 1970.
(13) Shaffer, M. J. and Dutt, G. R., Nitrification in Soil-Water
Systems: A Computerized Activated Complex Model. Proceedings
of the First World Congress on Water Resources. Chicago,
Illinois, 4:284-294, September 24-28, 1973.
(14) Shaffer, M. J., Dutt, G. R., and Moore, W. J., Predicting Changes
in Nitrogenous Compounds in Soil-Water Systems, Collected Papers
Regarding Nitrates in Agricultural Waste Water. U.S. Government
Printing Office, Washington, D.C., 1969. p. 15-28.
(15) Stockton, J. C. and Warrick, A. W., Spatial Variability of
Unsaturated Hydraulic Conductivity. SSSAP, 35:847-848, 1971.
(16) Dutt, G. R., Effect of Small Amounts of Gypsum in Soils on the
Solutes of Effluents. SSSAP. 28:754-757. 1964.
(17) Dutt, G. R., Prediction of the Concentration of Solutes in Soil
Solutions for Soil Systems Containing Gypsum and Exchangeable Ca
and Mg. SSSAP. 26:341-343. 1962.
(18) Dutt, G. R., and Anderson, W. D., Effect of Ca-Saturated Soils
on the Conductance and Activity of Cl , 50. = , and Ca'1"''. Soil
Science. 98: 377-382. 1964.
(19) Dutt, G. R. and Donneen, L. D., Predicting the Solute Composition
of the Saturation Extract from Soil Undergoing Salinization.
SSSAP. 27: 627-630. 1963.
(20) Dutt, G. R. and Tanji, K. K., Predicting Concentrations of Solutes
in Water Percolated Through a Column of Soil. Jr. Geophys. Res.
69:3437-3493. 1962.
(21) Dutt, G. R., Terkeltoub, R. W., and Ranschkolb, R. S., Prediction
of Gypsum and Leaching Requirements for Sodium Affected Soils.
Soil Science. 114:93-103. 1972.
226
-------
(22) Hanks, R. J. and Bowers, S. A., Numerical Solution of the Moisture
Flow Equation for Infiltration Into Layered Soils SSSAP
26:530-534. 1962.
(23) Campbell, G. S. "A Simple Method for Determining Unsatarated
Conductivity from Moisture Retention Data," Soil Science
Vol. 117, No. 6, pp. 311-314. 1974.
(24) Hillel, D. "Soil and Water, Physical Principles and Processes "
Academic Press, New York, 1971, pp. 288. '
227
-------
APPENDIX
The source code and sample output for the detailed model described
in this volume is available on microfische at cost from the authors,
Your request should be directed to:
Dr. Marvin J. Shaffer
Soil Scientist
United States Department of Interior
Bureau of Reclamation
Engineering and Research Center
P. 0. Box 25007
Denver Federal Center, Building 67
Denver, Colorado 80225
228
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
REPORT NO.
EPA-600/2-77-179e
2.
3. RECIPIENT'S ACCESSIOf*NO,
4. TITLE AND SUBTITLE
5. REPORT DATE
PREDICTION OF MINERAL QUALITY OF IRRIGATION RETURN
FLOW, VOLUME V, Detailed Return Flow Salinity and
and) Nutrient Simulation Model
August 1977 issuing date
6. PERFORMING ORGANIZATION CODE
7."AUTHOR(S)
Marvin J. Shaffer, Richard W. Ribbens
Charles W. Huntley
8. PERFORMING ORGANIZATION REPORT NO.
9. PERFORMING ORGANIZATION NAME AND ADDRESS
Bureau of Reclamation
Engineering and Research Center
Denver, Colorado 80225
10. PROGRAM ELEMENT NO.
1HB617
11. CONTRACT/GRANT NO.
EPA-IAG-D4-0371
12. SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Lab.
Office of Research and Development
U.S. Environmental Protection Agency
Ada, Oklahoma 74820
Ada, OK
13. TYPE OF REPORT AND PERIOD COVERED
Final
14. SPONSORING AGENCY CODE
EPA/600/15
15. SUPPLEMENTARY NOTES
VOLUMES I, II, III, IV
(EPA-600/2-77-179a thru 179d)
16. ABSTRACT
The development and evaluation of modeling capability to simulate and predict the
effects of irrigation on the quality of return flows are documented in the five
volumes of this report. The report contains two different modeling packages which
represent different levels of detail and sophistication. Volumes I, II, and IV
pertain to the model package given in Volume III. Volume V contains the more
sophisticated model. User's manuals are included in Volumes III and V.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.lDENTIFIERS/OPEN ENDED TERMS
COS AT I Field/Group
Mathematical Model, digital simulation,
scheduling, Irrigated land, Evapotrans-
piration, Agriculture, Agronomy, Water
pollution, Water loss.
Irrigation Return Flow
02 C/D
3. DISTRIBUTION STATEMENT
RELEASE TO PUBLIC
Unclassified
243
20. SECURITY CLASS (This page)
Unclassified
22. PRICE
EPA Form 2220-1 (9-73)
229
A U.S. GOVERNMENT PRINTING OFFICE: 1977- 757-056/6550
------- |