United States
Environmental Protection
Agency
Roberts Kerr Environmental Research      i-144
Laboratory
Ada OK 74820
Research and Development
Identification  and
Initial  Evaluation
of Irrigation
Return Flow
Models

-------
                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 poUution. 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-78-144
                                               July 1978
    IDENTIFICATION AND INITIAL EVALUATION OF
          IRRIGATION RETURN FLOW MODELS
                       by

                 Wynn R.  Walker
           Irrigation Hydrology Company
                  P.O. Box 1544
          Fort Collins, Colorado 80522
                Grant No. R804515
                 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 contents 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.

-------
                                  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.

     EPA1s 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 wastewaters with soil and
other natural systems; (c) develop and demonstrate pollution control tech-
nologies for irrigation return flows;  (d) develop and demonstrate pollution
control technologies for animal production wastes;  (e) develop and demon-
strate 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 pollu-
tion 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

-------
                                  ABSTRACT
     A broad based literature review was undertaken to identify studies that
had yielded digital computer models applicable to irrigation return flow
(IRF) systems.  The programs not listed in technical reports or papers were
requested from the various authors.  The results of this work are 43 computer
models applicable all or in part to the analysis of IRF's and their quality.
A brief evaluation of each model is given.

     IRF modeling technology is well developed theoretically but not
completely verified due to the large scale of the irrigation system.  Most
models remain in the research sphere and need to be redefined for the wider
utilization of planners.  Field data are generally not available to satisfy
the input requirements of most IRF models.  Accuracies of predictions need to
be determined against standardized conditions in order to further model
development and parameter sensitivities should be investigated to isolate the
most important field data.

     This report was submitted in fulfillment of Grant No. R-804515 to
Irrigation Hydrology Company, under the sponsorship of the U.S. Environmental
Protection Agency.  This report covers the period May 1, 1976 to November 30,
1977, and was completed as of November 30, 1977.

-------
                                  CONTENTS
                                                                      PAGE
Foreword ..............................
Abstract ..............................   iv
Figures ...............................   vi
Tables ...............................   yi
Acknowledgements ..........................  yii

1.   Introduction. ..... ....................    1
2.   Conclusions ..........................    3
3.   Recommendations ........................    5
4.   Methodology and Results ....................    7
          Model Identification ...................    "7
          Model Classification ...................    8
          Model Evaluation  .....................   12
          Assessment of IRF Modeling  ................   17

References  .............................   23
Bibliography  ............................   29
Appendices

     A.   Description of  Infiltration- Redistribution Models. ....   37

     B.   Description of  Unsaturated  Soil Solute Chemistry and
          Transport Models  ........  .  ............   51

     C.   Description of  Unsaturated  Zone Hydrology Models .....   62

     D.   Description of  Unsaturated  Zone Hydro-Quality Models  ...   75

     E.   Description of  Watershed Hydrology Models .........   103

     F.   Description of  Watershed Hydro-Quality Models .......   107

-------
NUMBER

  1.

  2.

  3.
                     LIST OF FIGURES



Schematic of the unsaturated zone in an irrigated system.

View of a large scale irrigation return flow system .  .  .

An illustration of a watershed which includes an
irrigated component 	
PAGE

 11

 13


 14
NUMBER
  1.
  2.
  3.
  4.
                               LIST OF TABLES
Classification of Irrigation Return Flow Models at
the Process Level 	
Classification of Irrigation Return Flow Models at
the Subsystem Level 	
Classification of Irrigation Return Flow Models at
the System Level	
                                                                      PAGE
                                                                       10
Summary of IRF Model Characteristics	16
                                      VI

-------
                              ACKNOWLEDGEMENTS

     The writer wishes to acknowledge the assistance of Mr. Stephen W. Smith,
President of Aqua Engineering Company and colleague at Colorado State
University for identifying the various computer programs reported in this
work.  Typing and drafting services were supplied by the Engineering Research
Center at Colorado State University.

     The writer is also indebted to Dr. Arthur G. Hornsby, Project Officer,
for cooperation during the course of the project and providing valuable
direction concerning its scope.
                                    VI1

-------
                                  SECTION 1

                                INTRODUCTION

     The use of computers in the study of the complex irrigation system
processes has provided a means for their rapid and intensive evaluation.  The
basis for the computer codes (or models) is the expression of the real system
as concise mathematical relationships.  The transformation between the
physical and mathematical forms has been achieved to various degress of
success depending on the expertise of the modeler, the objectives and funding
of the study, the required time and spatial scales of the investigation, and
the availability of data.  The major thrust of most models dealing with
irrigated agriculture has been associated with the excess water being applied
to croplands.  The volume of irrigation water needed to promote high crop
yields is almost universally less than the volumes diverted for these purposes.
Studies have been made to accurately determine the crop evapotranspiration so
that irrigation systems could be properly designed and operated.  The excess
moisture moving through the root zone has been investigated in many respects,
some of which include salt and fertilizer leaching and drainage requirements.
The quality of the irrigation wastewater  (called irrigation return flows,
IKP) has been of primary concern in recent years as efforts have been made to
maintain suitable water quality in all water resources.

     The principal interest of this work is the digital computer models which
deal with the environmental impact of irrigation return flows.  These are the
models which simulate the occurrence and magnitude of the flows, their
quality characteristics, and their effects on receiving waters.  Recognizing
that the IRF system is an aggregate of several complex hydrologic, chemical,
and biological processes, this group of models also includes simulations of
the basic processes.  The objective of this project was to identify and
initially evaluate available IRF computer programs.  The evaluation is
intended to be somewhat general, including a brief description of the model
scope, input requirements, characteristics of the computer code, expected
time and spatial resolutions, basic mathematical approach, and if possible a
note concerning previous applications, calibration, or verification.

     There is an underlying intent that these IRP models become more widely
utilized.  A substantial duplication of effort is evident in the formulation
and testing phases of various models, but this may be a minor problem com-
pared with the dissemination of the models to what might be called the
"applied users"  (planners, consultants, government agency personnel, etc.).
Most IRF models originate as research tools used to investigate the basic
behavior of irrigated agriculture.  Research interests seek an understanding
of the various hydro-quality mechanisms and their interactions enroute to
identifying the effectiveness of alternative management measures.  The

-------
research modeler attempts to simulate as much real system detail as is
possible.  Output is generally related to the time and spatial distribution
of parameters, and is generally compared to specially collected field data to
validate the approach and its assumptions.  Both time and spatial resolutions
are small  (many models study processes in a single vertical profile in terms
of minutes and extending no more than the period between irrigations).

     In contrast, there is a large number of people not so much interested in
understanding the problem as they are in getting it alleviated.  These are
the personnel charged with providing support to decision makers and adminis-
trators who recommend remedies or implementation policies to funding agencies
and legislatures.  Their view of the problem tends to be very macroscopic,
i.e., problems with inadequate field data for finely tuned research models, a
system of thousands or millions of hectares and time frames of months and
years.

     This study was organized to make an initial assessment of the existing
IRF modeling capability with a view towards improving the availability of the
programs for use by others.  In order that the reader might better understand
the selection of models, as will be discussed in later sections, several
points concerning the scope of the study should be made.  The models evaluated
in this report were identified from an intensive review of the technical
literature.  A large number of "models" are entirely mathematical or con-
ceptual, and are therefore not included.  Some models appear in the litera-
ture under different authors and applied to different situations.  For this
study, only the most recent or more complete version is included.  There
appear to be both studies in which a computer code is developed for general
use and those in which a program was written only for the purpose of the
specific project.  These latter models are not generally available for wide-
spread use because documentation has been omitted.  These programs which
would be difficult to use simply because no documentation is available or
because of rough coding are not part of this study.  It should also be noted
that not all models are actually available in the literature and many must be
requested from the authors.

     This report details some 43 individual models although the actual number
represented is much higher.  In addition to models having more than one
version or revision, some models are combinations of more specific models or
subroutines that may, in fact, often be utilized separately.  In this work,
the most integrated form of the models is presented.

     There are several major groundwater model collection and evaluation
efforts underway which are similar to this investigation.  To avoid dupli-
cating these efforts, groundwater models as such are not included herein.

-------
                                  SECTION 2

                                 CONCLUSIONS

     1.   The irrigation return flow system and many of its basic elements
are well simulated with the existing body of computer models.  However, the
number of individual models and their differences suggest that modeling
advances are not rapidly disseminated, particularly among related fields of
interest.  Except for some large scale watershed or river basin simulations,
few of the IRF models require large core and execution time commitments.

     2.   All of the programs evaluated were coded in either FORTRAN or CSMP.
Both languages have advantages to their users:  time and core requirements
are smaller for FORTRAN programs; programming is simpler and faster with
CSMP.  Many computers, particularly those other than IBM machines, will not
accommodate CSMP thereby limiting the widespread use of these programs.
FORTRAN versions are universally applicable.

     3.   Examination of the computer codes reveals a wide range of skills in
manipulating data, ordering calculations, and preparing output.  Models using
tape or disc drives are more difficult to apply on other computers because of
more variability in the format of these operations.  Card input is widely
used and offers little or no restrictions.  Few programs present output data
in an orderly and clear way.  In fact, many simply dump numbers without any
labeling.

     4.   The principal element in modeling any IRF system is the region from
the crop canopy to the water table.  This subgroup of IRF models is the best
developed and tested.  Unfortunately, available field data are insufficient
for most of these models so they are primarily research tools.  One of the
most severe data deficiencies is the description of the soil hydraulic prop-
erties, i.e., hydraulic conductivity as a function of moisture content, and
moisture content as a function of soil moisture potential.  If these data are
known, however, soil moisture flux is readily simulated.  The nitrogen cycle
remains the most elusive quality component of the unsaturated zone hydrology
primarily because it is an organic, biological series of transformations.
Evapotranspiration is a well modeled process when estimates are based on
climatological data.  Separation of transpiration and evaporation as well as
predicting the distribution of root water extraction is less developed.

     5.   Other subsets of the IRF models examined are also primarily slanted
toward research rather than planning applications.  Most if not all also
require at least some input data not generally available.  Consequently,
there is not an IRF modeling capability for evaluation of large scale
irrigated agriculture which can be utilized by the various nonresearch

-------
professionals.  This is a serious problem because it prevents an accurate
assessment of the effectiveness of alternative measures for managing the
quality of irrigation return flows.

     6.   Most models are written to be more general than actually required
in their initial use.  Consequently, there are often completely untested
loops within these programs.  Another user with a different set of conditions
may find problems in these heretofore unused loops.  There needs to be at
least two or more sets of data to serve as "standards" for the testing,
evaluation, and comparison of IRF models.

     7.   Although the models discussed in this report focus on the subject
of IRF quality, they are also valuable tools in other investigations.  The
basic processes simulated also provide the means for evaluating such topics
as energy conservation and storage in agricultural systems, maximizing yields
under drought conditions, disposal and impacts of organic wastes in agricul-
tural soils, and the feasibility of alternative irrigation systems.

-------
                                  SECTION 3

                               RECOMMENDATIONS

     1.   A library of the models identified in this study and others as they
become available should be organized at a central resource point for a wider
dissemination.  The library should include operational programs, a detailed
documentation for each model, and an example calculation,  provisions should
be made for periodic updating of the models and the acquisition of newly
developed programs.

     2.   There are many duplicative models.  Thus, each phase of the IRF
modeling technology should be condensed.  The relevant expertise should be
brought together with the purpose of developing a "second generation" series
of models by taking the best parts from all available models and producing a
series of models incorporating the strongest features of the existing
technology.

     3.   The CSMP models should be rewritten in FORTRAN for wider application
on other computer systems.

     4.   The research modelers need to determine the impact of variability
in their models.  Sensitivities to input  data, assumptions, and individual
parameters  should be  evaluated to delineate the kinds of information having
the highest priority  and  the limitations  of some mathematical approaches.
These  analyses  need to be made in a comparative format.  The sensitivity
studies need  to occur within the same  set of conditions  so that effect of
different modeling approaches can be ascertained.

     5.   Recommendation  4 requires two or more detailed data sets.  Because
of the high cost of field data collection, the available models need to be
subdivided  and  assessed one group at a time.  The first priority should be
given  to the  unsaturated  zone models.  Some data describing this system may
be already  available.  Consideration might also be given to a few completely
synthesized data sets representing conditions often encountered in field
situations.

     6.   This  largely research-oriented, modeling technology should be
redesigned  for  the much larger need involved in large scale planning studies.
Detailed models will  probably not be applicable, but should be used to test
the larger  scope models for their accuracy, dependability, and limitations.
Results of  this kind  of investigation  should be utilized to aid basic data
collection  efforts by detailing the most  important data  requirements and the
priority among  types  of data.

-------
     7.   Although there remain a number of modeling deficiencies, none
represent serious problems so far as irrigation return flow quality is
concerned.  A large enough basic research program is continuously operating
to remedy these deficiencies in a short time if reporting the respective
models can be encouraged.  The IRF models need to be visible enough to be
utilized in these new research activities so they will be subject to improve-
ment and expansion.

-------
                                  SECTION 4

                           METHODOLOGY AND RESULTS
MODEL IDENTIFICATION

     This project has taken a broad view of what constitutes an "irrigation
return flow model."  Actually, many of the computer programs were developed
for evaluation of individual or combined hydro-quality processes in irrigated
agriculture for entirely different objectives, but they are applicable to
studies which address the problems of IRF's.  Consequently, the models
assembled for this report are collectively denoted as IRF models because they
are usable and effective.  Another study concerned with crop production may
just as well refer to some of these same models as "production models"
because the same information is necessary.

     In recent years, several literature abstracting efforts have been made
to identify the pertinent irrigation return flow quality literature (Utah
State University Foundation  (1969) and Skogerboe et al. (1972, 1973, 1974,
1976, 1977, 1978)).  These references served as the initial starting point
for this study.  Key word indices were searched for such topics as computers,
computer models, models, numerical solutions, etc., which would delineate
reports where computer models had been developed or utilized.  Each potential
entry was obtained and evaluated in detail.  References cited in these
reports which indicated by their titles the further possibilities of model
development were also collected and reviewed.  This process was repeated
until no new literature on the subject of irrigation system models was
discovered.  Even then, personal contacts with colleagues revealed additional
modeling studies.  In addition, some literature entirely outside of agricul-
ture was identified in which models applicable to this subject were being
generated and used.  The final list of references dealing with the fairly
narrow topic of modeling totaled more than 200 citations.

     The principal interest  of this project was digital computer program or
code models that might be useful in investigating IRF systems.  Models which
are conceptual, completely mathematical, physical, or analog are omitted.
Many models are listed in reports and papers for the benefit of the potential
user and many are left to those sufficiently interested to request copies.
There is also a large body of models which are unavailable because they were
written for a specific purpose, not generalized to handle alternative forms
of the problem, and not documented or prepared for outside use.  Many of
these fortunately simulate what will be described later as process level
systems and are included in  other models.

-------
     After sifting the literature, the models requiring individual requests
were sought.  Letters were sent or telephone contacts were made with the
individual modelers asking for computer listings and available documentation.
Originally, an evaluation questionnaire was planned so that individual
modelers could describe their models for this work, but existing Federal
regulations prohibited this approach without special approval.  Response to
the written requests was generally poor.

     In earlier work, Walker  (1976) concluded that a great deal of model
duplication had occurred and that more effort to publicize various modeling
efforts would reduce this inefficient use of time and funding.  Later in this
section this conclusion will be supported.  However, in reviewing the litera-
ture there are several models which have been used by numerous other
researchers.  Some models have undergone several revisions and updates, or
have been slightly modified for another application.  The models referenced
in this report are generally the latest version, most complete, most applied,
or the ones most available.

     The models included herein do not represent an exhaustive listing since
there are undoubtedly models which this effort simply did not identify.
Nevertheless, the 43 models presented in the appendices encompass much of the
technology and certainly many of the best and most verified models available.

MODEL CLASSIFICATION

     In summarizing the models assembled for this study it is necessary to
segregate the programs into groups having similar scope and purpose.  To this
end, three major groups of two subclasses each were delineated:

     (1)  process level models
           (a)  infiltration-redistribution
           (b)  unsaturated soil solute chemistry and transport

     (2)  subsystem level models
           (a)  unsaturated zone hydrology
           (b)  unsaturated zone hydro-quality
     (3)  system level models
           (a)  watershed hydrology
           (b)  watershed hydro-quality

Following the subclass designations, the models are described in Appendices
A, B, C, D, E, and F.  A summary of the various models by major group is
given in Tables 1, 2, and 3.  In order to present the models more concisely
and to illustrate them more nearly as one would expect to use them in
modeling IRF systems, the models have been presented in their most aggregated
form.  The model by Shaffer, et al. (1977), for example, is presented as a
single system model even though it consists of seven individual subsystem
level models very much usable by themselves.
                                      8

-------
          TABLE 1.  CLASSIFICATION  OF  IRRIGATION RETURN
                    FLOW MODELS AT  THE PROCESS  LEVEL
      REFERENCE
INFILTRATION-
REDISTRIBUTION
 UNSATURATED  SOIL
 SOLUTE CHEMISTRY
   AND TRANSPORT

Bhuiyan, et al. (1971)
Brandt, et al. (1971)
Hillel (1977)
Lomen and Warrick (1974)
Van Der Ploeg and Benecke (1974)
Warrick (1974)
Warrick and Lomen (1976)
Wind and Van Doorne (1975)
Dugard and Reeves (1976)
Jurinak, et al. (1973)
Oster and McNeal (1971)
Rai and Franklin (1973)
Saxton, et al. (1977)
Tanji, et al. (1972)
(Appendix A)
A-l
A- 2
A- 3
A-4
A- 5
A- 6
A-7, A-8
A- 9






(Appendix B)








B-l
B-2.
B-3
B-4
B-5
B-6
           TABLE 2.  CLASSIFICATION OF IRRIGATION RETURN
                     FLOW MODELS AT THE SUBSYSTEM LEVEL
       REFERENCE
 UNSATURATED ZONE
    HYDROLOGY
   (Appendix C)
UNSATURATED ZONE
 HYDRO-QUALITY
 (Appendix D)
Feddes and Zaraday  (1977)
Gupta, et al.  (1977)
Hillel (1977)
Saxton, et al.  (1974)
Trava  (1975)
Neuman, et al.  (1975)

Beek and Frissel  (1973)
Davidson, et al.  (1975)
Davidson, et al.  (1977)
DeWit and Van  Keulen  (1975)
Frissel and Reiniger  (1974)
Hagin and Amberger  (1974)
Margheim  (1967)
Melamed, et al.  (1977)
Sroajstrla, et  al.  (1974)
       C-l
       C-2
    C-3,  C-4
       C-5
       C-6
       C-7
     D-10
                        D-l
                     D-2,  D-3
                       D-4
                        D-5
                   D-6, D-7, D-8
                        D-9
                        D-ll
                        D-12
                        D-l 3

-------
                TABLE 3.  CLASSIFICATION OF IRRIGATION RETURN
                          FLOW MODELS AT THE SYSTEM LEVEL
           P.EFEPENCE                   WATERSHED HYDROLOGY    WATERSHED
                                                            HYDRO-QUALITY
    	(Appendix E)	(Appendix F)

     Hillel  (1977)                             E-l
     Makkink and Van Heemst (1975)             E-2
Shaffer, et al. (1977)
Bureau of Reclamation (1977)
Crawford and Donigian (1973)
Dixon, et al. (1970)
Hill, et al. (1973)
Walker (1970)
F-l
F-2
F-3
F-4
F-5
F-6

Process Level Models

     The principal hydrologic processes in an irrigation return flow system,
illustrated in Figure 1, include transpiration, evaporation, infiltration-
redistribution, and root extraction.  Associated with each of these processes
are a number of chemical and biological transformations affecting the dis-
solved and adsorbed constituents in the soil.  The programs selected for this
report which might be classed as process level models are of two types:
(1) infiltration-redistribution; and  (2) unsaturated soil solute chemistry
and transport.  Other contributions to irrigation return flows which might be
categorized as process level topics, namely irrigation uniformity, conveyance
seepage, and field tailwater, are generally not modeled independently.
Irrigation uniformity has been recognized recently as an important input to
IRF simulations and several models that are available were not supplied to
the writer.  At least one uniformity model is being revised and validated at
Colorado State University, and should be available in the spring of 1978.
Conveyance seepage and field tailwater are most often addressed simply as
fractions of diverted or applied water.

     Root extraction, transpiration, and evaporation do not appear separately
in most modeling efforts.  Combined transpiration-evaporation (evapotranspi-
ration) can be determined from climatic conditions or evaluation of the
porous media physics.  Analyses based on climatic data are very simple and
commonly used so the relevant models are seldom presented as such.  Investi-
gations using soil physics approaches involve more than one process and are,
therefore, described in the subsystem level of modeling.

     The segregation of models into process, subsystem, and system levels may
reflect more of the writer's bias towards the modeling concept than is
justified.  However, process level models by themselves have only marginal
utility in evaluating IRF systems unless the modeler proposes to fabricate
new subsystem or system models.  Two important exceptions are to be noted.
First, deep percolation represents a major component of return flow.


                                     10

-------
                                        Transpiration
•
1 I

v_ v.
/
Root Extraction !
Upward Flow f
i » ' '-
V
1 Deep
r Percolation

nr Tnhle
*• ^v
-------
Analysis of deep percolation or leaching depends on either determining a root
zone mass balance or simulating the infiltration-redistribution process.  The
first alternative requires more data, is more applicable to large scale
situations, and is often the most desirable.  Evaluating infiltration-
redistribution involves shorter time intervals, less data (although more
sophisticated), and yields a closer approximation to the actual event.  The
second exception is the analysis of the chemical-biological aspects of soils
and soil solutes.  Moisture flow is often better approximated by assuming
steady-state conditions or that the transient nature may be represented by a
series of steady-state conditions.  In these cases, the quality aspects of
modeling are somewhat independent of the moisture phase and can be accurately
simulated by limiting the scope of the analyses.  These types of models are
also of benefit when computer time and core capacities are limited.  The
output of a moisture flux simulation can be input to chemical simulators
thereby significantly reducing the dimensionality of the program.

Subsystem Level Models

     In this report, the subsystem level of modeling is defined as the
combination of two or more process level simulations.  The main subsystem for
IRF studies is the field hydrology within the irrigated system itself.  In a
major watershed or river basin, there are other logical subsystems describing
precipitation-runoff relations, reservoir operations, in-stream processes
like dissolved oxygen behavior, groundwater interflows, etc.  Within the
irrigation system itself there may be a number of similar subsystems
reflecting different soils, crops, and irrigation management practices.

     As with any modeling classification, subsystem models vary in the scope
and dimension of their mathematics.  Some models are detailed treatments of a
small part of the irrigated agriculture subsystem shown in Figure 1, others
are more macroscopic, large scale and time resolved treatments of the Figure
2 representation.  The more detailed models involved basic physics and
chemistry whereas the macroscopic approaches are generally mass conservation
descriptions.

System Level Models

     As the scope and scale of IRF models increase, the variety of models
also becomes wider.  At the system level, there are two broad classes of
models included in this report.  The first is the models of the agricultural
watershed, shown in Figure 2, in which one or more unsaturated zone sub-
systems are linked with a groundwater or drainage subsystem simulation to
predict total IRF volumes and their time distributions.  The second set of
system models are those simulating river basin or subbasin hydrologies as
shown in Figure 3.  In these models, irrigated agriculture represents a small
fraction of the actual land area and the emphasis is on evaluating the impact
of irrigation on the system rather than quantifying its magnitude.

MODEL EVALUATION

     One of the purposes of this project was to make an initial assessment of
the various computer programs in order to identify the current modeling

                                      12

-------
            XG&ttgZ3r3^    L^VrX^^ %*\





•'.'.•.•••.V-'.'--  :->:-"-1:  ;•-; v.    v  :'•• -."  • v.
Root Zone
               ^_*J.  Wjter'_' Table
        Figure 2.   View of a large scale  irrigation return flow system.

-------
Figure 3.  An illustration of a watershed which includes an
           irrigated component (after Eakin et al. (1976)).

-------
capability.   Similar studies were examined to determine the alternatives for
describing the models.  After several revisions, the outline shown in the
appendices was selected.  The model descriptions are divided into the
following eight segments:

     (1)   Model grouping;
     (2)   Descriptive or developmental references;
     (3)   Scope of the model;
     (4)   Input requirements;
     (5)   Spatial and time scales;
     (6)   Structure of the code;
     (7)   Basic mathematical or analytical approach; and
     (8)   General comments.

     In describing various models, as much as possible of the information
presented originates in the cited references in order that descriptive errors
are minimized.  At the outset, it was desirable to maintain the length of the
summaries to one or two pages where possible.  This means that the phrases
used in describing various aspects of a model must alert the reader's own
understanding rather than trying to be exhaustive in detail.  Those unfamiliar
with a modeling area are encouraged to obtain and review the reference
material for a more in-depth assessment.  Much of the science behind these
models is beyond the training of the writer.  As a result, a description's
adequacy is highly dependent on the code documentation available with each
program.  A significant number of the models had little or no documentation
which could be used in this report in which case the writer attempted to
extrapolate from his own understanding.  Under all of these circumstances, it
is hoped that the essence of each model has been fairly represented.  A
summary of the various models' characteristics is given in Table 4.

Model Scope

     Under the heading of model scope, a few sentences were developed to
describe what the model  simulates and its expected output, including the use
of input dumps, error flags, and intermediate calculation results to make the
model more useable.  Most models print input data to insure that it has been
input properly.  Error flags and results of intermediate computations are not
generally included.

Input Requirements

     A list of control parameters and basic input data is included to
identify to the potential user the information required to utilize the models
effectively.  This list  is not exhaustive, but represents input data of most
importance.  Occassionally, a series of data is described in more general
terms or classified in order to reduce the length of the descriptions.

Spatial and Time Scales

     It is important  that the "size" of the modeled portion be identified in
comparison to the problem's real dimensions.  The interval or time scale is
also useful in selecting a model to evaluate an IRF problem since detail of

                                     15

-------
                                 TABLE 4.   SUMMARY OF IRF MODEL CHARACTERISTICS
cr>

MODEL
A-l
A- 2
A- 3
A-4
A- 5
A-6
A- 7
A-8
A- 9
B-l
B-2
B-3
B-4
B-5
B-6
C-l
C-2
C-3
C-4
C-5
C-6
C-7
D-l
D-2
D-3
D-4
D-5
D-6
D-7
D-8
D-9
D-10
D-ll
D-12
D-13
E-l
E-2
F-l
F-2
F-3
F-4
F-5
F-6

FLOWS EVALUATED
WATER SALINITY NITROGEN PESTICIDES
X
X
X
X
X
X
X
X
X






X
X
X (heat)
X
X
X
X
X
X
X
X
X (heat)
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X









X
X
X
X
X
X







X
X
X
X
X
X
X
X
X
X
X
X
X


X X
X X
X
X
X
X
STEADY (S)
OR
TRANSIENT (T)
T
T
T
T
T
T
T
T
T
T
S
S
S
T
S
T
T
T
T
T
S
T
T
T
S
T
T
T
T
T
T
T
S
T
T
T
S
T
T
T
S
S
S
FORTRAN
OR
CSMP
CSMP
Fortran
CSMP
Fortran
CSMP
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
CSMP
CSMP
Fortran
Fortran
Fortran
CSMP
Fortran
Fortran
Fortran
CSMP
CSMP
CSMP
CSMP
CSMP
CSMP
Fortran
Fortran
Fortran
CSMP
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
Fortran
RELATIVE
SIZE
Small
Medium
Small
Small
Medium
Small
Small
Small
Small
Medium
Small
Medium
Medium
Medium
Medium
Medium
Large
Small
Small
Large
Large
Medium
Medium
Large
Medium
Large
Medium
Medium
Medium
Medium
Large
Small
Medium
Medium
Medium
Medium
Medium
Large
Large
Large
Large
Large
Medium
VERIFICATION
REPORTED
Yes
Yes
No
Yes
Yes
No
No
No
No
No
Yes
Yes
Yes
Yes
Yes
Yes
Yes
No
No
Yes
Yes
Yes
No
Yes
Yes
Yes
No
No
No
No
No
No
No
Yes
Yes
No
No
Yes
Yes
Yes
Yes
Yes
Yes
DIMENSIONS
1
1, 2, 3
1
2
1, 2, 3
2
3
2
1
2
1
1
1
1
1
1
1
1
1
1
1
2
1
1
1
1
1
1
1
1
1
1
2
1
1
2
2
2
2
2
2
2
2

-------
the study is directly related to time resolution.  The dimension of the
simulation may also be provided in this section, thereby indicating the
assumptions made about the system's structure.

Computer Code Structure

     In these sections, the programming language, core requirements, expected
execution times, internal program linkage, and application comments are made.
Actually, it was discovered late in the study that core requirements and
execution times for the same program will vary substantially from computer to
computer.  Computer words have different byte lengths and execution times
depending on whether or not program compilation is necessary.  Consequently,
the writer tried to assess these parameters as they would be defined on the
CDC 6400 computer at Colorado State University.  A major problem occurs when
the code is written in a CSMP format because much of the program then becomes
part of the basic computer software and storage-time requirements are not
known.

Basic Mathematical Approach

     It is difficult to express in one or two paragraphs what a complex
simulation model accomplished and in what manner.  Key phrases and words are
used as much as possible to convey to the reader a great deal more than is
actually written, and if possible these phrases were taken from the documen-
tation of the code or reference material.  The assumptions in a model are at
least partially apparent to the already informed reader from the description
of what is calculated.  Often classical equations in the literature are named
without reference.  These are noted in the text and generally identify the
analytical approach taken.

General

     Of general interest to modelers is where the model can be obtained,
verification that has occurred, the stage of development and future applica-
tions, and special problems that might be encountered in its use.  For the
most part, this section contains comments on verification and acquisition.

ASSESSMENT OF IRF MODELING

General

     If irrigation return flows are to be minimized in an area, the unused
diversion should approximate the water necessary to maintain a favorable salt
balance in the cropland root zone.  The capacity of conveyance and drainage
systems should be reduced to the peak seven to ten day average consumptive
use rate and the minimal leaching fraction, respectively.  Large capital
investments in remodeling the irrigation and drainage facilities will have
to be made in advance and each irrigator will require substantial retraining.
(It can also be successfully argued that the agricultural scientists and
engineers would need significantly better practical or "hands-on" training in
order to assist the irrigator.)  The quality of the flows proceeding to the
sequential downstream users will still be somewhat degraded because of the

                                     17

-------
evaporative water depletion, but the effect will be minimized.  This objective
may appear to be difficult to achieve, however, in at least one instance in
western Colorado, there are already authorized plans to invest $100-$200
million into a salinity control project involving 25,000 hectares with the
goal that salt loading might be reduced by something on the order of 400,000
metric tons annually.  In this particular case and many others of a similar
nature, years of investigation preceded the decision for implementation.  The
dimension of the problem was defined and the effectiveness of the alternative
improvements were predicted.  The IRF models were instrumental in accom-
plishing this work and they should continue to be the most important tools
available to planners and researchers as they seek to remedy other irrigation
return flow problems.

     Modeling is the art of model development and testing, but the term
"model" itself covers a broad range.  In this report, the group of models of
interest is the so called "computer models," or those whose mathematical basis
is executed by digital computers.  These models are developed from concise
mathematical expressions of the real system being evaluated, based on the
best available evidence of how the system behaves under existing conditions.
If the mathematics are correct and general, the same model will predict or
simulate the behavior of the real system under changed conditions.  Thus, it
is possible to expand the information on the system to areas where data have
not been collected, thereby yielding a greater understanding of the problems
involved in its management.

     IRF models can be described in general terms by noting some of the
classical methods of model description.  None of the models collected for this
work or noted in the literature are stochastic in nature.  In other words, IRF
models are made deterministic by assuming the parameters and variables are
uniquely related and are not distributed on the basis of probability.  The
types of variables and parameters in IRF models are both lumped and distri-
buted depending primarily on the modeler's assumptions regarding homogeneity.
Process and many subsystem models have distributed parameter characteristics
to separately describe the variation in hydrologic conditions.  Some of the
subsystem and most of the system level models have lumped parameters because
input is usually the average conditions.  The solution techniques applied to
IRF model mathematics are almost invariably numerical because of the inherent
complexity being evaluated.  Analytical solutions are possible after compara-
tively significant simplifying assumptions are made, and these solutions are
often utilized to test a more general numerical model.  The continuity or time
resolution of IRF models also varies with the scale of the simulation.
Process models generally approach a continuous simulation by utilizing time
increments as small as one-tenth day whereas the system models are generally
more discrete as their time increments are on the order of months.  Time
resolutions are also affected by the dynamic versus steady-state approach
being applied.  Small time events such as soil moisture flux require a dynamic
or time-dependent simulation to evaluate changes.  When the only answer sought
is the total or aggregate flux from the root zone, a steady-state analysis is
often used.  Because the steady-state system is generally easier to simulate
some models actually "simulate" the dynamic or transient characteristics with
a series of steady-state solutions  (the CSMP approaches are among the most
refined of these approaches).

                                     18

-------
     The IRF system is very complex and not always adequately defined.
Studies range in the scope in which the IRF system must be investigated.
Consequently, there are numerous models covering the same simulation scope
and reflect the understanding and interpretation of the system by the modeler
himself.  What is not well known at this time is how much difference can be
expected when applying these models to the same conditions and which models
to apply under various conditions.

Computer Coding

     There are basic differences among various modelers in the computer
coding approach taken; however, the most significant differences are related
to the scope of the model.  Process and small subsystem models tend to be
straightforward and condition-specific while system level models are usually
written to be much more generally applicable.  As a result, a great deal more
input is required in system level models to define the problem structure and
the alternative means for its evaluation.

     It appears to be standard practice to utilize subroutines to separate
calculations in most FORTRAN models.  Data are about equally transferred by
the call statement and through common dimensioning.  These practices have an
advantage that should not be overlooked: the abstraction and individual use
of important subroutines by other modelers.  Another advantage in actual
coding is the relative ease of debugging.  A disadvantage is that programs
tend to be longer and require more core and time resources.

     The use of the IBM simulation language, CSMP, is said to simplify the
coding process, but reported comparisons suggest execution times may be
greater.  The CSMP models cannot be applied to other computers as FORTRAN
models can, and this limitation should be strongly considered.

     Documentation within the computer code is very useful to users and
should be considered.  Most IRF models are poorly documented within the code
itself.

     Probably the biggest problem in the computer codes themselves is that
they contain "subtle" programming errors, a misplaced parenthesis, mislabeled
statements, redundant looping, etc.  This is particularly a problem with
large programs.  When a model is written and tested, it is unlikely that the
modeler's own problem will require looping through each alternative route in
the program.  The output is corrected for the gross mistakes because they are
apparent; the subtle mistakes are much less easily corrected.  Thus, when
another user applies the model to another set of conditions, difficulty is
often encountered.  It would be very useful to have several standard data
sets on which various computer programs could be tested to evaluate their
accuracy, codes, and parameter sensitivities.

The State-of-the-Art

     The reader can review the appendix material and several recent reports
to decide for himself what the level of technology or state-of-the-art is in
IRF modeling.  In fact, most reports and papers review the pertinent literature

                                      19

-------
in preparation to presenting their information.  Consequently, these next few
paragraphs will address the IRF modeling state-of-the-art from a different
perspective.

     When one contemplates the body of IRF models available today and
considers the approach that will be taken in reducing the environmental
impact from these flows, it becomes apparent that:   (1) the existing models
are almost entirely research oriented rather than applied or planning tools;
and (2) the most important future need is not in research application but in
planning and implementing water quality management strategies.  The models
require data that are simply not available outside a research investigation
and are not easily or inexpensively collected.  Areawide planners are not
using these models in developing strategies for pollution abatement and risk
the possibility that these plans will not be cost-effective.  There needs to
be an effort made to bridge this gap as soon as possible.  The existing
models need to be thoroughly evaluated again, their parameter sensitivities
quantified, and alternative simplifying assumptions assessed.  A new modeling
technology designed for planning and large scale evaluation efforts needs to
be generated from the research model along with a quantitative description of
their functions and expectations.  The scientific principles surrounding the
IRF system should be reviewed by the consultants and government agencies
involved in water quality planning.

     Many IRF models demonstrate a noticeable lack of interdisciplinary
involvement and yet very good application of interdisciplinary principles.
This indicates that literature reviews and professional associations are
broad based and active.  The use of available computer programs noted by the
literature, however, is not generally practiced; rather, most modelers
utilize the concepts and mathematics but not the codes themselves.  The
detailed soil chemistry simulation reported by Dutt, et al.  (1972) and Shaffer,
et al. (1977) is possibly the best exception.  This model has been utilized
in studies in several universities and federal research groups.  Nevertheless,
the central tendency has been to develop individual programs from the common
mathematical core.  After having reviewed the codes of some 43 of these
models, it is the writer's attitude that receding has been the easiest route
because of poor program documentation.

     As the reader will become aware in reviewing the appendix model
descriptions, there are some modeling phases better developed than others.
Estimating the potential and actual evapotranspiration at the soil or crop
surface from climatological data is both well developed and tested, but many
IRF models employ energy balance or combination equations which require
substantially more data than are generally available.  The use of the equa-
tions relying on available data will not result in significant loss of
accuracy if applied properly (although the minimal time increments will be
extended).   Evaluation of the unsaturated zone hydrology often requires that
transpiration and evaporation be separated.  Approaches using leaf area index
or soil physics have been developed and appear in some of the models included
herein.  These simulations are more complicated and thus less accurate than
climatic functions.  The major difficulty lies in root extraction and growth
processes.   Most models assume uptake is linearly related to'the root mass in
a soil layer and that the root mass distribution is defined apriori.  The

                                      20

-------
more complete models attempt to include submodels of root growth and activity
and then introduce stress effects as the limitations of hydraulic conductivity
near the roots.  In these cases, the summation of uptake is equated to actual
transpiration as a check on mass balance and as a means to iteratively distri-
bute the uptake.  An interesting sidelight here is the differences in which
stress effects are handled in the two approaches,  With the climatic estimates,
stress is an empirical relationship based on total soil moisture depletion
whereas the soil physics approach actually limits transpiration by treating
it as a sink term in the moisture flow equation.  With adequate data, it
should be possible to numerically improve the climatic approach stress
function by linking the two methods together.

     Within the unsaturated soil profile, infiltration and redistribution of
an irrigation or rainfall event occurs over a relatively short period of time
(2 to 3 days in many soils).  In these situations, evapotranspiration is a
small fraction of the moisture movements and is generally neglected.  It
should be noted that the advent of high frequency irrigation systems remove
the acceptability of these assumptions.  Nearly every model evaluated in this
study encompassing unsaturated flow uses the same form of the one phase
moisture flow  equation derived by combining the time-depth dependent forms
of the Darcy equation and conservation of mass.  The soil moisture character-
istics  (hydraulic conductivity versus soil moisture content versus soil
matrix potential) are required for all of these models.  These are data
universally unavailable outside of study areas.  If a rapid in-situ field
technique were available, these models would add a great deal to the capa-
bility of consultants and planners.  The unsaturated soil moisture models
solve the nonlinear partial differential equations  (for the transient cases)
with finite difference, finite element, or the  CSMP simulation of sequential
steady-state approximations.  Numerical dispersion is corrected in a number
of ways.  Thus, infiltration-redistribution is  comparatively well modeled
with the most  severe limitation being available  data.

     If the model scope is expanded to where the unsaturated zone subsystem
is viewed in the context of its chemical constituent, the reliability is
somewhat diminished.  The mass transport of the  salts, including diffusion
and dispersion, is comparatively well simulated  by miscible displacement
theory.  Even  the exchange reactions with the  soil and solid salt phase, the
solution-precipitation of slightly soluble salts, and the formation of ion
pairs in the soil solution have been modeled successfully.  Some soils
contain very large volumes of solid phase gypsum and lime salts adding to the
complexity of  the salinity modeling, but again verified modeling efforts are
reported in the literature.  The same is not completely true of solutes
involving nitrogen compounds and certain volatile pesticides.  Nitrogen
transformations are generally simulated with first-order kinetics using rate
coefficients reported in the literature or developed in small scale labora-
tory studies.   Denitrification continues to be  the most difficult segment of
the nitrogen system to model while ammonium and  nitrate transport are the
easiest.  Pesticides often adsorb to the soil matrix and volatilize during
their tenure in the soil profile, both processes are evaluated in the
existing models.
                                      21

-------
     When the scale of IRF modeling increases beyond the unsaturated zones of
the irrigated field hydrology, a great deal of change occurs in the modeling
scope.  Most of these models are steady-state, lumped parameter types.  Their
scale represents the simulation of entire fields, valleys, and river basins
and the answers sought are those reflecting how irrigation return flows
affect the quality of a major flow resource.  Input-output or mass balance is
the predominant driving force in these system level models.  Thus, the
detailed simulations at the subsystem and process level are often not
directly applicable to large scale models.  The problem of how to achieve the
proper vertical integration in IRF modeling remains unsolved.  There are
models of major hydrologic subbasins that include a somewhat detailed sub-
program of soil physics and soil chemistry, but add very little, if any,
additional modeling capability at this scale because of the difficulty in
defining averages on so fine a scale.  In a large irrigated region, the
"average" is a concept rather than a reality to be used as data for simu-
lating the area.  Many IRF models do not encounter this disparity in resolu-
tion but it is worthwhile to note nonetheless.

     A final point should be made concerning the state-of-the-art in IRF
modeling.  Verification has been a major effort in the development of process
and subsystem models.  The system scale models are generally no more than
calibrated.  The difference being that model output when compared to measured
data  (verification) is possible when the scope is small, whereas the adjust-
ment of parameters until the predicted output is the same as the measured
data  (calibration) is the extent of what can be done in simulating large
scale systems.
                                     22

-------
                                 REFERENCES

Beek, J. and M. J. Frissel.  1973.  Simulation of Nitrogen Behavior in Soils.
     Centre for Agricultural Publishing and Documentation, Wageningen, The
     Netherlands.  67 pp.

Bhuiyan, S. I., E. A. Hiler, C. H. M. Van Bavel, and A. R. Aston.  1971.
     Dynamic Simulation of Vertical  Infiltration Into Unsaturated Soils.
     Water Resources Institute, Texas ASM University, College Station,
     Texas.  23 pp.

Brandt, A., E. Bresler, N. Diner, I. Ben-Asher, J. Heller, and D. Goldberg.
     1971.  Infiltration from a Trickle Source:  I.  Mathematical Models.
     Soil Sci. Soc. Amer. Proc.  35:675-682.

Bresler, E.  1973.  Simultaneous Transport of Solute and Water Under
     Transient, Unsaturated Flow Conditions.  Water Resources Research
     9(9):975-986.

Bresler, E., J. Heller, N. Diner, J. Ben-Asher, A. Brandt, and D. Goldberg.
     1971.  Infiltration from a Trickle Source:  II Experimental Data and
     Theoretical Productions.  Soil Sci. Soc. Amer. Proc., 35:683-689.

Bureau of Reclamation.  1977.  Prediction of Mineral Quality of Irrigation
     Return Flow, Volume III, Simulation Model of Conjunctive Use and Water
     Quality for a River System or Basin.  User's Manual.  EPA-600/2-77-179c.
     Robert S. Kerr Environmental Research Laboratory, Office of Research and
     Development, U.S. Environmental Protection Agency, Ada, Oklahoma.
     August.  285 pp.

Childs, S. W., and R. S. Hanks.  1975.  Model of Soil Salinity Effects on
     Crop Growth.  Soil Sci. Soc. Amer. Proc.  39:617-622.

Crawford, N. H., and A. S. Donigian.  1973.  Pesticide Transport and Runoff
     Model for Agricultural Lands.  EPA-660/2-74-013, Office of Research and
     Development, U.S. Environmental Protection Agency, Washington, D.C.
     December.  317 pp.

Davidson, J. M., G. A. Brusewitz, D. R. Baker and A. L. Wood. 1975.  Use of
     Soil Parameters for Describing Pesticide Movement Through Soil.
     EPA-660/2-75-009.  National Environmental Research Center, U.S.
     Environmental Protection Agency, Corvallis, Oregon.  May.  149 pp.
                                      23

-------
Davidson, J. M., D. A. Graetz, P. S. C. Rao, and H. M. Selim.  1977.
     Simulation of Nitrogen Movement, Transformations, and Uptake in the
     Plant Root Zone.  Draft Report to the U.S. Environmental Protection
     Agency, Office of Research and Development, Athens, Georgia.  139 pp.

DeWit, C. T., and H. Van Keulen.  1975.  Simulation of Transport Processes in
     Soils.  Centre for Agricultural Publishing and Documentation.  Wageningen,
     The Netherlands.  100 pp.

Dixon, N. P., D. W. Hendricks, A. L. Huber, and J. M. Bagley.  1970.
     Developing a Hydro-quality Simulation Model.  Report PRWG 67-1.  Utah
     Water Research Laboratory, Logan, Utah.  June.  193 pp.

Donigian, A. S., and N. H. Crawford.  1976.  Modeling Pesticide and Nutrients
     on Agricultural Lands.  EPA-600/2-76-043.  Office of Research and
     Development, Environmental Research Laboratory, U.S. Environmental
     Protection Agency, Athens, Georgia.  February.  211 pp.

Dugard, J. D., and M. Reeves.  1976.  Material Transport Through Porous
     Media:  A Finite-element Galerkin Model.  Environmental Sciences
     Division Publication 733, Oak Ridge National Laboratory, Oak Ridge,
     Tennessee.  March.  201 pp.

Dutt, S. R., M. S. Shaffer, and W. S. Moore.  1972.  Computer Simulation
     Model of Dynamic Bio-physicochemical Processes in Soils.  Technical
     Bulletin 196, Agricultural Experiment Station, University of Arizona,
     Tucson, Arizona.  October.  101 pp.

Eakin, T. E., D. Price, and Harrill.  1976.  Summary Appraisals of the Nation's
     Groundwater Resources—Great Basin Region.  Geological Survey Professional
     Paper 813-G.  U.S. Government Printing Office, Washington, D.C.  37 pp.

Feddes, R. A., and H. Zaraday.  1977.  Numerical Model for Transient Water
     Flow in Nonhomogeneous Soil-root Systems with Groundwater Influence.
     Proceedings of IFIP Working Conference on Modeling and Simulation of
     Land, Air, and Water Resource Systems, Ghent.  September.  18 pp.

Frissel, M. J., and P. Reiniger.  1974.  Simulation of Accumulation and
     Leaching in Soils.  Centre for Agricultural Publishing and Documenta-
     tion, Wageningen, The Netherlands.  116 pp.

Gupta, S. K., K. K. Tanji, D. R. Nielsen, J. W. Biggar, S. C. Simmons, and
     J. L. Maclntyre.  1977.  Field Simulation of Soil-water Movement with
     Water Extractions.  Water Science and Engineering Paper 4013.  University
     of California, Davis, California.  May.  95 pp.

Hagin, J., and A. Amberger.  1974.  Contribution of Fertilizers and Manures
     to the N- and P-Load of Waters, A Computer Simulation Model.  Report to
     the Deutsche Forschung Gemeinschaff.  January.  123 pp.
                                     24

-------
Hill,  R.  W. ,  E. K. Israelsen, and J. P. Riley.  1973.  Computer Simulation of
     the Hydrologic and Salinity Systems Within the Bear River Basin.
     PRWG 104-1.  Utah Water Research Laboratory/ Utah State University,
     Logan, Utah.  122 pp.

Hillel,  D.  1977.  Computer Simulation of Soil-Water Dynamics.  International
     Development Research Centre, Ottawa, Canada.  214 pp.

Huber, A. L., E. K. Israelsen, R. W. Hill, and J. P. Riley.  1976.  Basin
     Simulation Assessment Model Documentation and Users Manual.  PRWG 201-1.
     Utah Water Research Laboratory, Utah State University, Logan, Utah.
     August.   30 pp.

Jensen,  M. E., C. E. Franzoy, and D. C. N. Robb.  1970.  Scheduling
     Irrigations Using Climate-crop-soil Data.  J. Irr. and Drain. Div.,
     ASCE, 96:25-38.

Jurinak, J. J., S. Lai, and J. J. Hassett.  1973.  Cation Transport in Soils
     and Factors Affecting Soil Carbonate Solubility.  EPA-R2-73-235,  Office
     of Research and Monitoring, U.S. Environmental Protection Agency,
     Washington, D.C.  May.  87 pp.

Kincaid, D. C., and D. F. Heerman.  1974.  Scheduling Irrigations Using a
     Programmable Calculator.  ARS-NC-12.  Agricultural Research Service,
     U.S. Department of Agriculture.  February.  55 pp.

Lomen, D. 0.  and A. W. Warrick.  1974.  Time-dependent Linearized
     Infiltration:  II.  Line Sources.  Soil  Sci. Soc. of Amer. Proc.,
     38:568-572.

Margheim, G.  A.  1967.  Predicting the Quality of Irrigation Return Flows.
     Unpublished M.S. Thesis, Civil Engineering Department, Colorado State
     University, Fort Collins, Colorado.  December.  57 pp.

Makkink, G. F., and H. D. J. Van Heemst.  1975.  Simulation of the Water
     Balance of Arable Land and Pastures.  Centre for Agricultural Publishing
     and Documentation, Wageningen, The Netherlands.  79 pp.

Melamed, J. D., R. J. Hanks, and L. S. Willardson.  1977.  Model of Salt Flow
     in Soil with a Source-sink Term.  Soil Sci. Soc. Amer. Proc., 41:29-33.

Narasimhan, V. A., and E. K. Israelsen.  1975.  A Water-land Use Management
     Model for the Sevier River Basin, Phases I and II.  Information Bulletin
     24.  Utah Water Research Laboratory, Utah State University, Logan, Utah.
     September.  44 pp.

Neuman, S. P., R. A. Feddes, and E. Bresler.   1975.  Finite Element Analysis
     of Two-dimensional Flow in Soils Considering Water Uptake by Roots:  I.
     Theory,  II.  Field Applications, Soil Sci. Soc. Amer. Proc., 39:224-237.
                                      25

-------
Nimah, M. D. , and P., .T. Hanks.  1973.  Model for Estimating Soil Water, Plant,
     and Atmospheric Interrelations:  I. Description and Sensitivity, and II.
     Field Test of Model.  Soil Sci. Soc. Amer. Proc., 37:522-527 and 37:528-532.

^toucj., u. LI. ~i^i ts. .LJ. Hi-weal.  j.9"/j..  Computation of Soil Solution
     Composition Variation with Water Content for Desaturated Soils.  Soil
     Sex. £oc. .-iiner. Proc., 35:^36-442.

Oster, J. D. and J. D. Rhoades.  1975.  Calculated Drainage Water Composition
     and Salt Burden Resulting From Irrigation with River Waters in the
     Western United States.  J. Environ. Qual., 4:73-79.

Ostrowski, J. T.  1976.  Behavior of Groundwater Subject to Irrigation of
     Effluent—A Case Study.  M.S. Thesis, Civil Engineering Department,
     Maryland University, College Park, Maryland.  121 pp.

Rai, D., and W. T. Franklin.  1973.  Program for Computing Equilibrium
     Solution Composition in CaCO  and CaSO  Systems from Irrigation Water
     Compositions.  Water Management Technical Report No. 29.  Colorado State
     University, Fort Collins, Colorado.  October.  42 pp.

Saxton, K. E., H. P. Johnson, and R. H. Shaw.  1974.  Modeling Evapotranspira-
     tion and Soil Moisture.  Trans. ASAE, 17:673-677.

Saxton, K. E., G. E. Schuman, and R. E..Burwell.  1977.  Modeling Nitrate.
     Movement and Dissipation in Fertilized Soils.  Soil Sci. Soc. Amer.
     Proc., 41:265-271.

Shaffer, M. S., R. W. Ribbens, and C. W. Huntley.  1977.  Prediction of
     Mineral Quality of Irrigation Return Flow, Volume V.  Detailed Return
     Flow Salinity and Nutrient Simulation Model.  EPA-600/2-77-179e.  Robert
     S. Kerr Environmental Research Laboratory, Office of Research and
     Development, U.S. Environmental Protection Agency, Ada, Oklahoma.
     August.  228 pp.

Skogerboe, G. V., V. T. Sahni, and W. R. Walker.  1972.  Selected Irrigation
     Return Flow Quality Abstracts 1968-1969.  EPA-R2-72-094.  Office of
     Research and Monitoring, U.S. Environmental Protection Agency,
     Washington, D.C.  October.  211 pp.

Skogerboe, G. V., W. R. Walker, D. J. Meyer, R. S. Bennett.  1973.  Selected
     Irrigation Return Flow Quality Abstracts 1970-1971.  EPA-R2-73-271.
     Office of Research and Monitoring, U.S. Environmental Protection Agency,
     Washington, D.C.  June.  285 pp.

Skogerboe, G. V., W. R. Walker, R. S. Bennett, and B. J. Zakely.  1974.
     Selected Irrigation Return Flow Quality Abstracts 1972-1973.  EPA-660/2-
     74-049.  Office of Research and Development, U.S. Environmental Protec-
     tion Agency, Washington, D.C.  June.  409 pp.
                                      26

-------
Skogerboe, G. V., W. R. Walker, and S. W. Smith.  1976.  Selected Irrigation
     Return Flow Quality Abstracts, 1974.  EPA-600/2-76-019.  Office of
     Research and Development, U.S. Environmental Protection Agency,
     Washington, D.C. March.  220 pp.

Skogerboe, G. V., S. W. Smith, and W, R. Walker.  1977.  Selected Irrigation
     Return Flow Quality Abstracts, 1975.  EPA-600/2-77-094.  Robert S. Kerr
     Environmental Research Laboratory, Office of Research and Development,
     U.S. Environmental Protection Agency, Ada, Oklahoma.  May.  249 pp.

Skogerboe, G. V., S. W. Smith, and W. R. Walker.  1978 (In Press).  Selected
     Irrigation Return Flow Quality Abstracts, 1976.  Robert S. Kerr Environ-
     mental Research Laboratory, Office of Research and Development, U.S.
     Environmental Protection Agency, Ada, Oklahoma.  310 pp.

Smajstrla, A. G., D. L. Reddell, and E. A. Hiler.  1974.  Simulation of
     Miscible Displacement in Soils.  Paper 74-2015 presented at the Annual
     Meeting of ASAE, Stillwater, Oklahoma.  June.  31 pp.

Tanji, K. K., L. D. Doneen, G. V. Ferry, and R. S. Ayars.  1972.  Computer
     Simulation Analysis on Reclamation of Salt-affected Soils in San Joaquin
     Valley, California.  Soil Sci. Soc. Amer. Proc., 36:127-133.

Trava, T. L.  1975.  Allocating Available Irrigation Water on the Farm.
     Unpublished Ph.D. Thesis, Department of Agricultural Engineering,
     Colorado State University, Fort Collins, Colorado.  December.  107 pp.

Utah State University Foundation.  1969.  Characteristics and Pollution
     Problems of Irrigation Return Flow.  Robert S. Kerr Environmental
     Research Laboratory, Office of Research and Development, U.S. Environ-
     mental Protection Agency, Ada, Oklahoma.  May.  237 pp.

Van Der Ploeg, R. R., and P. Benecke.  1974.  Unsteady, Unsaturated,
     n-Dimensional Moisture Flow in Soil.  Soil Sci. Soc. Amer. Proc.,
     38:881-885.

Walker, W. R.  1970.  Hydrosalinity Model of the Grand Valley.  Unpublished
     M.S. Thesis, Civil Engineering Department, Colorado State University,
     Fort Collins, Colorado.  August.  94 pp.

Walker, W. R.  1976.  Assessment of Irrigation Return Flow Models.  EPA-600/2-
     76-0219.  Robert S. Kerr Environmental Research Laboratory, Office of
     Research and Development, U.S. Environmental Protection Agency, Ada,
     Oklahoma.  October.  75 pp.

Warrick, A. W.  1974.  Time-dependent Linearized Infiltration.  I. Point
     Sources.  Soil Sci. Soc. Amer. Proc., 38:383-386.

Warrick, A. W., and D. 0. Lomen.  1976.  Time-dependent Linearized
     Infiltration:  III.  Stip and Disc Sources.  Soil Sci. Soc. Amer. Proc.,
     40:639-643.
                                      27

-------
Wind, G. P., and W. Van Doorne.  1975.  A Numerical Model for the Simulation
     of Unsaturated Vertical Flow of Moisture in Soils.  J. of Hydrology,
     24:1-20.
                                      28

-------
                                BIBLIOGRAPHY

Abendt, R. W.  1975.  Models in Water Quality Planning.  Ecological Modelling,
     1(3):205-217.

Bailey, G. W., R. R. Swank, Jr., and H. P. Nicholson.  1974.  Predicting
     Pesticide Runoff from Agricultural Land:  A Conceptual Model.  Journal
     of Environmental Quality, 3(2):95-102.

Bassett, D. L.  1972.  Mathematical Model of Water Advance in Border
     Irrigation.  TRANSACTIONS of the ASAE, 15(5):992-995.

Bliesner, R. D., R. J. Hanks, L. R. King, and L. S. Willardson.  1976.
     Effects of Irrigation Management on the Quality of Irrigation Return
     Flow.  Soil Sci. Soc. Amer. Jour., 41(2):424-428.

Boast, C. W.  1973.  Modeling the Movement of Chemicals in Soils by Water.
     Soil Sci., 115:224-230.

Bouwer, H.  1969.   Infiltration of Water Into Nonuniform Soil.  ASCE,
     Proceedings, Journal of Irrigation and Drainage Division, Paper 6937,
     95(IR4):451-462.

Bouwer, H.  1976.   Infiltration Into Increasingly Permeable Soils.  ASCE,
     Journal of Irrigation and Drainage Division  (102):129-136.

Braester, C., G. Dagin,  S. Neuman, and D.  Zaslavsky.   1971.  A Survey of the
     Equations and  Solutions of Unsaturated Flow in Porous Media.  First
     Annual Report, Project A10-SWC77.  Hydraulic Eng. Lab., Technion, Haifa,
     Israel.  176 pp.

Branson, R. L., P.  F. Pratt, J. D. Rhoades, and J. D.  Oster.  1975.  Journal
     of Environmental Quality, 4(1):33-40.

Bresler, E.  1967.  A Model for Tracing Salt Distribution in the Soil Profile
     and Estimating the  Efficient Combination of Water Quality and Quantity
     Under Varying  Field Conditions.  Soil Sci., 104:227-233.

Bresler, E.  1972.  Interacting Diffuse Layers in Mixed Mono-Divalent Ionic
     Systems.  Soil Sci. Soc. Amer. Proc., 36(6):891-396.

Bresler, E.  1973.  Anion Exclusion and Coupling Effects in Nonsteady
     Transport Through Unsaturated Soils.  I.  Theory.  Soil Sci. Soc.
     Amer. Proc., 37:663-669.
                                      29

-------
Bresler, E., and R. J. Hanks.  1969.  Numerical Method for Estimating
     Simultaneous Flow of Water and Salt in Unsaturated Soils.  Soil Sci.
     Soc. Amer. Proc., 33:827-832.

Bresler, E., and D. Yaron.  1972.  Soil Water Regime in Economic Evaluation
     of Salinity in Irrigation.  Water Resources Research, 8(4):791-800.

Bresler, E., and A. Laufer.  1974.  Anion Exclusion and Coupling Effects in
     Nonsteady Transport Through Unsaturated Soils:  II.  Laboratory and
     Numerical Experiments.  Soil Sci. Soc. Amer. Proc., 38 (2):213-218.

Bresler, E., and D. Russo.  1975.  Two-dimensional Solute Transfer During
     Nonsteady Infiltration:  Laboratory Test of Mathematical Model.  Soil
     Sci. Soc. Amer. Proc., 39(3):585-587.

Brutsaert, W. F.  1971.  A Functional Iteration Technique for Solving the
     Richards Equation Applied to Two-dimensional Infiltration Problems.
     Water Resources Research, 7:1583-1596.

Brutsaert, W. F., and R. N. Weisman.  1970.  Comparison of Solutions of a
     Nonlinear Diffusion Equation.  Water Resources Research, 6(2):642-644.

Burman, R. D., and R. D. Black.  1970.  The Inference of Intake and
     Hydraulic Roughness Parameters from Plot Runoff Using Kinematic Wave
     Theory.  TRANSACTIONS of  the ASAE, .13(4):479-481.

Chen, L. H., B. K. Huang, and W. E. Splinter.  1969.  Developing a Physical-
     Chemical Model for a Plant Growth System.  TRANSACTIONS of the ASAE,
     12(5):698-702.

Cheng, R. T., and C. Y. Li.  1973.  On the Solution of Transient Free-surface
     Flow problems in Porous Media by the Finite Element Method.  Journal of
     Hydrology, 20(l):49-63.

Chidley, T. R. E., and J. G. Pike.  1970.  A Generalized Computer Program for
     the Solution of the Penman Equation for Evapotranspiration.  Journal of
     Hydrology, 10(1):75-89.

Childs, E. C., and M. Bybordi.  1969.  The Vertical Movement of Water in
     Stratified Porous Material.  1.  Infiltration.  Water Resources Research,
     5(2):446-459.

Childs, S. W., and R. J. Hanks.  1975.  Model to Predict the Effect of Soil
     Salinity on Crop Growth.  Soil Sci. Soc. Amer. Proc., 39:617-622.

Davidson, J. M., D. R. Baker, and G. H. Brusewitz.  1975.  Simultaneous
     Transport of Water and Adsorbed Solutes Through Soil Under Transient
     Flow Conditions.  TRANSACTIONS of the ASAE, 18(3):535-539.

Day, P. R., and J. N. Luthin.  1956.  A Numerical Solution of the Differential
     Equation of Flow for a Vertical Drainage Problem.  Soil Sci. Soc. Amer.
     Proc., 20:443-447.

                                      30

-------
Drake, R. L., and C. P. Peterson.  1971.  Application of a Local Similarity
     Concept in Solving the Vertical Subsurface Flow Problem.  Water Resources
     Research, 7(5):1241-1255.

Eagleman, J. R.  1971.  An Experimentally Derived Model for Actual
     Evapotranspiration.  Agricultural Meteorology, 8(4/5):385-394.

Enfield, C. G., and D. C. Shew.  1975.  Comparison of Two Predictive
     Nonequilibrium One-dimensional Models for Phosphorus Sorption and
     Movement Through Homogeneous Soils.  Journal of Environmental Quality,
     4(2):198-202.

Freeze, R. A.  1969.  The Mechanism of Natural Groundwater Recharge and
     Discharge-1.  One-dimensional, Vertical, Unsteady, Unsaturated Flow
     Above a Recharging or Discharging Groundwater Flow System.  Water
     Resources Research, 5(1):153-171.

Garder, A. D., D. W. Peaceman, and A. L. Pozzi.  1964.  Numerical Calculation
     of Multidimensional Miscible Displacement by the Method of Characteris-
     tics.  Soc. Petrol. Engr. J., 4:26-36.

Green, D. W., H. Dabiri, C. F. Weinaug, and R. Prill.  1970.  Numerical
     Modeling of Unsaturated  Groundwater Flow and Comparison of the Model to
     a Field Experiment.  Water  Resources Research, 6(3):862-874.

Griffin, R. A., and J. J. Jurinak.  1973.  Test of a New Model for the
     Kinetics of Adsorption-Desorption Processes.  Soil Sci. Soc. Amer.
     Proc., 37(6):869-872.

Gupta, S. C.  1972.  Model for Predicting Distribution of Salt and Water in
     Soils.  Ph.D. Dissertation, Utah State University, Logan, Utah.

Hanks, R. J., and S. A. Bowers.  1962.  Numerical Solution of the Moisture
     Flow Equation for Infiltration Into Layered Soils.  Soil Sci. Soc. Amer.
     Proc., 26:530-534.

Hanks, R. J., A. Klute, and E. Bresler.  1969.  A Numerical Method for
     Estimating Infiltration, Distribution, Drainage, and Evaporation of
     Water from Soil.  Water  Resources Research, 5:1064-1069.

Hart, W. E.  1972.  Subsurface Distribution of Nonuniformly Applied Surface
     Waters.  TRANSACTIONS of the ASAE, 15(4):656-661.

Hornberger, G. M., I. Remson, and A. A. Fungaroli.  1969.  Numeric Studies
     of a Composite Soil Moisture Groundwater System.  Water Resources
     Research, 5(4):797-802.

Hornberger, G. M., and I. Remson.  1970.  A Moving Boundary Model of a
     One-dimensional Saturated-Unsaturated, Transient Porous Flow System.
     Water Resources Research, 6(3):898-905.
                                     31

-------
Hornsby, A. G.  1973.  Prediction Modeling for Salinity Control in Irrigation
     Return Flows.  EPA-R2-73-168.  U.S. Environmental Protection Agency,
     Corvallis, Oregon.  55 pp.

Hunt, B. W.  1973.  Dispersion in Nonuniform Seepage.  Journal of the
     Hydraulics Division, ASCE Paper 9553, 99(HY2):295-299.

Huntley, C. W.  1974.  Abbreviated Users Manual for Program EVPCOM.  Division
     of Planning Coordination, Bureau of Reclamation, Engineering and
     Research Center, Denver, Colorado.  100 pp.

Huntley, C. W.  1974.  Abbreviated Users Manual for Program ETBYMO.  Division
     of Planning Coordination, Bureau of Reclamation, Engineering and
     Research Center, Denver, Colorado.  36 pp.

Hurley, P. A.  1968.  Predicting Return Flows from Irrigation.  Journal of
     Irrigation and Drainage Division, ASCE Paper 5838, 94(IRl):41-48.

James, L. G., and C. L. Larson.  1974.  Modeling Infiltration and
     Redistribution of Soil Water During Intermittent Applications.
     Presented at 1974 Winter meeting of the ASAE, December 10-13, Chicago,
     Illinois.  16 pp.

Jamieson, D. G., and C. R. Amerman.  1969.  Quick-return Subsurface Flow.
     Journal of Hydrology, 8(2):122-136.

Kastanek, F.   1971.  Numerical Simulation Technique for Vertical Drainage
     from a Soil Column.  Journal of Hydrology, 14 (3/4):213-232.

King, L. G., and R. J. Hanks.  1973.  Irrigation Management for Control of
     Quality of Irrigation Return Flow.  EPA-R2-73-265.  U.S. Environmental
     Protection Agency, Washington, D.  C.  306 pp.

Lai, S. H., and J. Jurinak.  1971.  Numerical Approximation of Cation
     Exchange in Miscible Displacement  Through Soil Columns.  Soil Sci. Soc.
     Amer. Proc., 35:894-899.

Lai, H., and J. Jurinak.  1972.  Cation Adsorption in One-dimensional Flow
     Through Soils:  A Numerical Solution.  Water Resources Research,
     8:99-107.

Lindstrom, F. T., R. Haque, and W. R. Coshow.  1970.  Adsorption from
     Solution:  III.  A New Model for the Kinetics of Adsorption-Desorption
     Processes.  J. Phy. Chem., 74:495.

Marino, M. A.  1974.  Numerical and Analytical Solutions in a Finite,
     Adsorbing Porous Medium.  Water Resources Bulletin, 10(1):81-90.

Marino, M. A.  1974.  Distribution of Contaminants in Porous Media Flow.
     Water Resources Research, 10(5):1013-1018.
                                     32

-------
Mehran, M., and K. K. Tanji.  1974.  Computer Modeling of Nitrogen
     Transformations in Soils.  Journal of Environmental Quality, 3(4):391-395.

Mein, R. G., and C. L. Larson.  1973.  Modeling Infiltration During a Steady
     Rain.  Water Resources Research, 9(2):384-394.

Misra, C., D. R. Nielsen, and J. W. Biggar.  1974.  Nitrogen Transformations
     in Soil During Leaching:  I.  Theoretical Considerations.  Soil Sci.
     Soc. Amer. Proc., 38:289-293.

Misra, C., D. R. Nielsen, and J. W. Biggar.  1974.  Nitrogen Transformations
     in Soil During Leaching:  II.  Steady State Nitrification and Nitrate
     Reduction.  Soil Sci. Soc. Amer. Proc., 38:294-299.

Molz, F. J.  1972.  Simulation of  Post-Irrigation Moisture Movement.  ASCE,
     Journal of Irrigation and Drainage Division, 98(IR4):523-532.

Molz, F. J.  1974.  Practical Simulation  Models of the Subsurface Hydrologic
     System with Example Applications.  Water Resources Research Institute,
     Auburn University, Auburn, Alabama,  Bulletin 19.  46 pp.

Molz, F. J., and I. Remson.  1970.  Extraction Term Models of Soil Moisture
     Use by Transpiring Plants.  Water Resources Research, 6:1346-1356.

Molz, F. J., and I. Remson.  1971.  Application of an Extraction Term Model to
     the Study of Moisture Flow to Plant  Roots.  Agron. J., 63:72-77.

Morey,  R. V., and J. R. Gilley.   1973.  A Simulation  Model for Evaluating
     Irrigation Management Practices.  TRANSACTIONS  of the ASAE, 16(5):
     979-983.

Neuman, S. P.  1973.   Saturated-Unsaturated  Seepage  by Finite Elements.
     ASCE, Journal  of  the Hydraulics  Division,  99(HY12):2233-2250.

Novak,  L.  T., and D. C. Adriano.   1975.   Phosphorus  Movement in  Soils:
     Soil-Orthophosphate Reaction Kinetics.  Journal of Environmental
     Quality, 4 (2):261-266.

Oster,  C., J. Sonnichsen, and  R.  Jaske.   1970.  Numerical Solution  to  the
     Convective Diffusion Equation.   Water Resources Research,6(6):1746-1752.

Passioura,  J. B., and  I. R.  Cowan.  1968. Solving the Nonlinear Diffusion
     Equation for the  Radial Flow of  Water to Roots.  Agricultural
     Meteorology,  5:129-134.

Perez,  A.  I., W.  C.  Huber,  J.  P.  Heaney,  and E. E.  Pyatt.  1972.  A Water
     Quality Model  for Conjunctive Surface-Groundwater System:   An  Overview.
     Water Resources Bulletin,  8(5):900-908.

Philip,  J.  R.   1955.   Numerical  Solution of Equations of the Diffusion Type
     with Diffusivity  Concentration Dependent.  Faraday  Soc.  Trans.,
      51:885-892.

                                      33

-------
Remson, I., A. A. Fungaroli, and G. M. Hornberger.  1967.  Numerical Analysis
     of Soil Moisture Systems.  Proc. Amer. Soc. Civil Engr., J. Irr. and
     Drain. Div., IRS:153-166.

Reuss, J. O., and J. M. Geist.  1970.  Prediction of the Nitrogen Requirements
     of Field Crops.  Part I:  Theoretical Models of Nitrogen Release.
     Agronomy Journal 62(3):381-384.

Routson, R. C., and R. J. Serne.  1972.  Experimental Support Studies for the
     PERCOL and Transport Models.  Battelle Pacific Northwest Laboratories
     Report BNWL-1719.  Richland, Washington.  71 pp.

Routson, R. C., and R. J. Serne.  1972.  One-dimensional Model of the
     Movement of Trace Radioactive Solutes Through Soil Columns:  The PERCOL
     Model.  BNWL-1718, Battelle Pacific Northwest Laboratories.  Richland,
     Washington.  December.

Rubin, J.  1968.  Theoretical Analysis of Two-dimensional, Transient Flow of
     Water in Unsaturated and Partly Unsaturated Soils.  Soil Sci. Soc. Amer.
     Proc., 32:602-615.

Seginer, I., and J. Morin.  1970.  A Model of Surface Crusting and Infiltration
     of Bare Soils.  Water Resources Research, 6(2):629-633.

Serne, R. J., and R. C. Routson.  1972.. Users' Manual - PERCOL:  A Soil-.
     Waste Chemical Equilibrium Model.  BNWL-1720.  Battelle Pacific North-
     west Laboratories.  Richland, Washington.  December.

Serne, R. J., R. C. Routson, and D. A. Cochran.  1973.  Experimental Method
     for Obtaining PERCOL Model Input and Verification Data.  Battelle
     Pacific Northwest Laboratories Report BNWL-1721.  Richland, Washington.
     38 pp.

Shah, D. B., G. A. Coulman, L. T. Novak, and B. G. Ellis.  1975.  A
     Mathematical Model for Phosphorus Movement in Soils.  Journal of
     Environmental Quality, 4(1):87-92.

Shastry, J. S., L. T. Fan, and L. E. Erickson.  1973.  Nonlinear Steady-State
     Model of Isotropic Fractionation Accompanying Nitrogen Transformations
     in Soil.  Soil Sci. Soc. Amer. Proc., 38(2):315-322.

Skaggs, R. W.  1975.  A Water Management Model for High Water Table Soils.
     Presented at the 1975 Winter Meeting of the ASAE, December 15-18,
     Chicago, Illinois.  27 pp.

Skopp, J., and A. W. Warrick.  1974.  A Two-Phase Model for the Miscible
     Displacement of Reactive Solutes in Soils.  Soil Sci. Soc. Amer. Proc.,
     38(4):545-550.

Smith, R. E.  1972.  Border Irrigation Advance and Ephemeral Flood Waves.
     ASCE, Journal of the Irrigation and Drainage Division, 98(IR2)289-307.
                                     34

-------
Stammers, W. N., O. C. Igwe, and H. R. Witeley.  1973.  Calculation of
     Evaporation from Measurements of Soil Water and the Soil Water
     Characteristic.  Canadian Agricultural Engineering, 15(1):2-5.

Stubb, W. J.  1966.  Infiltration and Redistribution of Water in Vertical
     Columns of Loam Soil.  Soil Sci. Soc. Amer. Proc., 30:553-558.

Tanji, K. K.  1970.  A Computer Analysis on Leaching of Boron From
     Stratified Soil Columns.  Soil Science, 110:44-51.

Tanji, K. K.  1976.  A Conceptual Hydrosalinity Model for Predicting Salt
     Load in Irrigation Return Flows.  Proceedings of the International
     Conference on Managing Saline Water for Irrigation:  Planning for the
     Future, Texas Tech University, Lubbock, Texas,  pp. 49-70.

Tanji, K. K., and L. O. Doneen.  1966.  A Computer Technique  for Prediction
     of CaCO  Precipitation-HCO~ Salt Solutions.  Soil Sci. Soc. Amer. Proc.,
     30:53-56.

Tanji, K. K., G. R. Dutt, J. L. Paul, and I. 0. Doneen.  1967.  Quality of
     Percolating Waters.  II.  A Computer Method for Predicting Salt
     Concentrations in Soils at Variable Moisture Contents.   Hilgardia,
     38:307-318.

Tanji, K. K.,  I. O. Doneen, and J. L. Paul.  1967.  Quality of Percolating
     Water.  III.  Predictions on  the Quality  of Water Percolating Through
     Stratified Substrata by Computer Analysis.  Hilgardia, 38:319-347.

Taylor,  L.  S.,  and J. N. Luthin.   1969.  Computer Methods for Transient
     Analyses  of Water Table Aquifers.  Water  Resources Research,  5:144-152.

Thomas,  A.  W.,  E. G.  Kruse, and H. R. Duke.  1974.  Steady Infiltration from
     Line Sources Buried in Soil.  TRANSACTIONS of the ASAE,  17(1):125-129.

Thomas,  J.  L.f  J. P.  Riley, and E. K. Israelsen.  1972.  A Hybrid  Computer
     Program  for Predicting the Chemical Quality of Irrigation Return Flows.
     Water  Resources  Bulletin, 8(5):922-934.

Todorovic,  P.   1970.  A Stochastic Model of  Longitudinal Diffusion in Porous
     Media.  Water  Resources Research,  6(1):211-222.

Todsen,  M.  1973.   Numerical Studies  of Two-dimensional Saturated-Unsaturated
     Drainage  Models.  Journal of  Hydrology, 20(4):311-326.

Van Genuchten,  M.  T., J. M. Davidson, and P. J. Wierenga.  1974.   An
     Evaluation of Kinetic  and Equilibrium Equations  for the  Prediction of
     Pesticide Movement Through Porous  Media.   Soil Sci. Soc. of Amer. Proc.,
      38(1): 29-35.

Van Keulen, H., and D. G. E. M. Van  Beek.  1971.  Water Movement in Layered
      Soils.   A Simulation Model.   Neth. J. Agr. Sci.,  19:138-153.
                                      35

-------
Walker, A.  1974.  A Simulation Model for Prediction of Herbicide Persistence.
     Journal of Environmental Quality, 3(4):396-401.

Walter, M. F., T. S. Steenhuis, G. D. Bubenzer, and J. C. Converse.  1974.
     Evaluation of a Soil Nitrate Transport Model.  Presented at 1974 Winter
     Meeting of the ASAE,  December 10-13, Chicago, Illinois.  18 pp.

Warrick, A. W.  1974.  Solution to the One-dimensional Linear Moisture Flow
     Equation with Water Extraction.  Soil Sci. Soc. Amer. Proc., 39(4):
     573-576.

Warrick, A. W., and D. Kirkham.  1969.  Two-dimensional Seepage of Ponded
     Water to Full Ditch Drains.  Water Resources Research, 5(3):685-693.

Warrick, A. W., J. W. Biggar, and D. R. Nielsen.  1971.  Simulation of Solute
     and Water Transfer for an Unsaturated Soil.  Water Resources Research,
     7:1216-1225.

Wei, C. and R. W. Jeppson.  1971.  Finite Difference Solutions of
     Axisymmetric Infiltration Through Partially Saturated Porous Media.
     Water Research Lab. Report PRWG 59C-6.  Utah State University, Logan,
     Utah.  April.  64 pp.

Whisler, F. O., and A. Klute.  1965.  The Numerical Analysis of Infiltration
     Considering Hysteresis into a Vertical Soil Column at Equilibrium Under
     Gravity.  Soil Sci. Soc. Proc., 29:489-494.

Whisler, F. D., K. K. Watson, and S. J. Perrens.  1972.  Numerical Analysis
     of Infiltration into Heterogeneous Porous Medium.  Soil Sci. Soc. Amer.
     Proc., 36:868-874.

Wind, G. P.  1972.  A Hydraulic Model for the Simulation of Nonhysteretic
     Vertical Unsaturated Flow of Moisture in Soils.  Journal of Hydrology,
     15 (3):227-246.
                                     36

-------
                                 APPENDIX A

              DESCRIPTION OF INFILTRATION-REDISTRIBUTION MODELS


MODEL:                                A-l

REFERENCES:
Bhuiyan, S. I., E. A. Hiler, C. H. M. Van Bavel, and A. R. Aston.  1971.
Dynamic simulation of vertical infiltration into unsaturated soils.  Water
Resources Institute, Texas ASM University, College Station, Texas.  23 pp.

MODEL SCOPE:
     Simulation is of the infiltration into and through a soil profile.
Output includes not only the time and spatial distributions of the infiltrated
water but the time distributed prediction of the intake relationship.

INPUT REQUIREMENTS;

     Input data include moisture content, hydraulic conductivity, and soil
suction relations, and initial moisture distribution.  Time increments are
calculated from user specified error tolerances.

SPATIAL AND TIME SCALES:
     Water movement in a  one-dimensional unsaturated soil profile is
analyzed.  Time periods for calculations are in seconds, while output
is aggregated up to hours.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The unsteady infiltration of water into a soil medium is analyzed by
assuming a ponded surface, uniform initial soil moisture distribution, and
an impermeable  layer at  the  lower soil boundary.  The unsaturated soil is
then divided into layers for the following two step iterative process.
First, the flux between  layers is determined using mass conservation, Darcy's
Law, and the hydraulic potential between layers.  Second, the incremental
water content of each layer  is then determined by integrating the net flux
                                      37

-------
rate with a fourth order Runge-Kutta method.  This updated moisture content
is then the initial condition for subsequent iteration.

GENERAL:
     The authors note that the program can be modified to handle buried water
sources, root extraction, and multi-layered soils.  CSMP code is listed in
the reference cited above.  Verification trials are reported.
                                      38

-------
MODEL:                                A-2

REFERENCES:

Brandt, A., E. Bresler, N. Diner, I. Ben-Asher, J. Heller, and D. Goldberg.
1971.  Infiltration from a trickle source:  I. Mathematical models.  Soil
Sci. Soc. Amer. Proc., 35:675-682.

Bresler, E.r J. Heller, N. Diner, I. Ben-Asher, A. Brandt, and D. Goldberg.
1971.  Infiltration from a trickle source:  II. Experimental data and
theoretical predictions.  Soil Sci. Soc. Amer. Proc., 35:683-689.

MODEL SCOPE:
     Simulation is of the moisture distribution pattern under a trickle
irrigation source.  Analysis allows for transient conditions and involves
both a plane flow and cylindrical flow description.  Hysteresis is neglected.
Output includes the time and spatial distribution of soil moisture,
evaporation, and deep percolation.

INPUT REQUIREMENTS;

     The input data include saturation moisture content, initial distribution
of soil moisture, diffusivity and hydraulic conductivity versus moisture
content, discharge of the trickle source and its duration and frequency.  A
number of control parameters set the model in either plane or cylindrical
flow coordinates, number of grid points, number and size of time increments,
convergence tolerances, number of irrigation cycles, and other geometrical
variables.

TIME AND SPATIAL SCALES:
     Analysis of the trickle  system is based on the application and
redistribution of the irrigation water.  Calculations are based on small
intervals  (minutes or seconds) but the basic interval is the irrigation
frequency.  Model is a  two-dimensional treatment of the root zone.

COMPUTER CODE STRUCTURE:
     The code is a Fortran  IV program using  several subroutines with call
statement data transfer  and several functions using common data linkage.
Core requirements and  execution time are not discussed.

BASIC MATHEMATICAL APPROACH:
     The transient infiltration and redistribution from a trickle irrigation
source is simulated with a  conventional diffusion type water flow equation
using a finite difference technique  (a combination of the noniterative
alternating directions  (ADI) procedure with a Newton-Raphson technique).
Boundary conditions describe the trickle emission and evaporation, and the
soil conditions  surrounding the expected wetting zone.
                                      39

-------
GENERAL;

     The model has been verified against field data and utilized by several
other users.  The code must be requested from the authors.
                                    40

-------
MODEL:                                A-3

REFERENCES:

Hillel, D.  1977.  Computer Simulation of Soil-Water Dynamics.  International
Development Research Centre, Ottawa, Canada,  pp. 35-60.

MODEL SCOPE:
     Analysis is of the isothermal evaporation from a bare soil water under
fluctuating evaporation potential.  Soil is homogeneous and isotropic.  Out-
put includes the soil surface evaporation and the time-spatial distribution
of soil moisture in the soil profile.

INPUT REQUIREMENTS;

     The model requires definition of the system geometry  (depth of soil
profile and number of layers, thickness of layers if not uniform), moisture
content, hydraulic conductivity, and moisture suction relations, initial
moisture distribution, time distributed evaporation demand at the surface,
and the air dry matric potential.  If hysteresis is to be considered, addi-
tional soil moisture-matric suction tables must be supplied to define the
sorption-desorption relationships.

TIME AND SPATIAL SCALES:
     The calculations occur on one second  intervals with aggregation to days
if desired.  The  spatial  scale is a one-dimensional treatment of a 1-2
meter unsaturated soil profile.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The soil profile  is  divided  into  compartments or layers.  Movement
between compartments follows  Darcy's Law  and  conservation of mass.  The
computation procedure  follows a sequential or iterative evaluation of com-
partment moisture content and then  calculation of flux between compartments
based on averaged gradients and hydraulic conductivity.  The evaporation rate
from the soil surface  is  set  equal  to  the potential rate when the top layer's
soil moisture level is above  the  air-dry  value.  If below, the evaporation
rate is set equal to the  flux into  the upper  layer.

GENERAL:
     The model was  numerically  tested with and without hysteresis, but no
 field verification  is  reported.  Results  show significant effects and might
 be usefully  applied elsewhere in evaluating hysteresis.  The CSMP code is
 listed  in  the reference.
                                      41

-------
MODEL;                                A-4

REFERENCES;

Lomen, D.O. and A.W. Warrick.  1974.  Time-dependent linearized infiltration:
II.  Line Sources.  Soil Sci. Soc. Amer. Proc., 38:568-572.

MODEL SCOPE:

     Analysis of infiltration and redistribution of irrigation water from a
line source such as multi-outlet trickle irrigation lateral is presented.
Moisture flux or matrix flux potential is the primary simulation variable.
Flux is expressed in relation to time and distance from the sources.

INPUT REQUIREMENTS;

     Input to the model includes the following:  (1) number of line sources;
(2) line at which computations start; (3) distance between lines;  (4) real
dimension grid width, depth, and time; and (5) soil hydraulic properties.
Some input is required for checking of various numerical functions.

SPATIAL AND TIME SCALES;

     Spatial scale varies from the width of an irrigated field to the
unsaturated depth.  Time scale encompasses the duration of an irrigation.

COMPUTER CODE STRUCTURE;

     Written in Fortran IV and using 26 k bytes of storage and approximately
1-2 minutes of central processor time, program is easily adaptable to various
computers.  A main program coordinates functions through call and common
information exchange.

BASIC MATHEMATICAL APPROACH;

     The computer code solves the line source infiltration problem on the
assumption of step input of irrigation water and a linearized relationship
between hydraulic conductivity and moisture content.  These assumptions allow
analytic solution to the matric flux potential relation as described in the
reference cited above.

GENERAL:

     Program copy requests should be addressed to the authors cited above.
Verification has been made with field data but not reported at the time of
this writing.
                                      42

-------
MODEL:                                A-5

REFERENCES;

Van Der Ploeg, R. R. and P. Benecke.  1974.  Unsteady, unsaturated,
n-dimensional moisture flow in soil.  Soil Sci. Soc. Amer. Proc., 38:881-885.

MODEL SCOPE:
     Simulation is of one, two, or three dimensional unsteady, unsaturated
flow in soils.  Output includes the time and spatial distribution of moisture
contents and fluxes.  Four individual models are actually available:
(1) 1-dimensional absorption;  (2) 3-dimensional absorption  (spherical);
(3) 2-dimensional infiltration from a line-source; and  (4)  3-dimensional
infiltration from a point source.  The analytical structure is the same in
every case so these models will be presented collectively.

INPUT REQUIREMENTS;

     The potential user must divide the soil profile into "compartments" with
a parameter.  The multi-dimensional analysis requires definition in each
dimension.  The hydraulic properties of the soil such as moisture content
dependent conductivity and diffusivity must be defined.  Initial conditions
through the profiles as well as the distributed input at the surface are also
supplied.

SPATIAL AND TIME SCALES:
     The analysis is carried out in  fractions of an hour with integrated
output in intervals of hours or days.  Spatial scales are limited to a small
region examined in one, two, or three dimensions.

COMPUTER CODE STRUCTURE:
     All four codes are written in the  IBM  360/365 version of CSMP.
Execution time for the simplest model is about 10 minutes.

BASIC MATHEMATICAL APPROACH:
     The compartment approach described herein is discussed for other CSMP
models included in this work.  The profile is divided into segments depending
on the dimensionality of  the model and defined with respect to their initial
hydraulic properties.  Then, the ratio of moisture movement between compart-
ments is calculated on the basis of Darcy's Law and conservation of mass by
assuming the hydraulic gradient, diffusivity, and hydraulic conductivity are
the average between compartments.  For small time increments, the moisture
content in each compartment as computed by summing the fluxes into and out of
the compartment.  The procedure is then repeated for each compartment until
the desired time of examination is reached.  These models utilize a fourth
order Runge-Kutta integration method to sum events and do not rely on
constant time or step sizes.


                                     43

-------
GENERAL;

     Models have been calibrated and verified against field data in several
places.  The code listing is contained in the reference.
                                      44

-------
MODEL;                                A-6

REFERENCES;

Warrick, A.W., 1974.  Time-dependent linearized infiltration. I. Point
Sources.  Soil Sci. Soc. Amer. Proc., 38:383-386.

MODEL SCOPE;

     Simulation is of infiltration and redistribution of irrigation water
from a point source such as a trickle irrigation emitter.  Input dump and
system test are provided.  Model will simulate single or multiple outlets.

INPUT REQUIREMENTS;

     The program first requires a series of arrays to test the numerical
integration and exponential functions.  The number of points, strength or
discharge, dimensionless depth, radius, and time points are defined.  Corre-
sponding to each point defined, the associated initial conditions are read.
Use of the program requires that both the hydraulic conductivity and the
moisture content versus pressure head relationships be known.

SPATIAL AND TIME SCALES;

     The analysis involves the two-dimensional simulation of soil moisture
flux into a soil of known properties during a single or multiple irrigation
event.  The linearized analysis allows multiple points to be evaluated.

COMPUTER CODE STRUCTURE;

     The code consists of a main program and functions served by common and
call data transfer.  Approximately 25 k bytes and 10-20 central processor
seconds are required to operate the model.

BASIC MATHEMATICAL APPROACH;

     Analytic solutions developed and reported by the author based on the
assumptions of step-input of irrigation water and linear relationship between
hydraulic conductivity and moisture content.  See model A-4.

GENERAL;

     This analysis is most applicable where moisture content changes are
small such as a high frequency or trickle irrigation systems.  Computer code
must be requested from the original source.  Limited verification has been
accomplished.
                                     45

-------
MODEL;                                A-7

REFERENCES:

Warrick, A.W., and D.O. Lomen.  1976.  Time-dependent linearized infiltration:
III.  Stip and disc sources.  Soil Sci. Soc. Amer. Proc., 40:639-643.

MODEL SCOPE:

     Simulation is of disk sources of moisture input at the soil surface.
Analysis is three-dimensional.  Parameters included in simulation are prima-
rily matric fluxes as a function of depth in soil and distance  (radius) from
the source origin.  Input dump and error flags are incorporated.

INPUT REQUIREMENTS;

     The user first specifies the number of values of disc radius to be
evaluated, and the number of dimensionless radii, depths, and times.  Dimen-
sionless parameters require knowledge of the hydraulic conductivity-pressure
head relationship and the hydraulic conductivity-moisture content function.
The values of radius, depth and time  (dimensionless) representing initial
soil moisture conditions are then read in at each grid point established.

SPATIAL AND TIME SCALE;

     The model is a three-dimensional simulation of soil moisture flux
during an irrigation interval.  System solution allows the additive feature
so that sequential and parallel disks can be evaluated.

COMPUTER CODE STRUCTURE;

     The program incorporates a main program and function routines using
common and call data transfer.  Code is written in Fortran IV and uses 40 k
bytes and 5-10 central processor seconds per analysis.  Program is usable in
most computer systems as currently written.

BASIC MATHEMATICAL APPROACH;

     The generally accepted moisture flow equation for the time-dependent
case becomes a linear differential equation if the hydraulic conductivity is
assumed linearly related to moisture content.  This assumption allows the
matric potential to be described by analytic solutions which are discussed in
the reference.  (See model A-4.)

GENERAL;

     Code must be requested from the authors.
                                     .46

-------
MODEL;                                A-8

REFERENCES:

Warrick, A.W., and D.W. Lomen.  1976.  Time-dependent linearized infiltration:
III.  Strip and disc sources.  Soil Sci. Soc. Amer. Proc., 40:639-643.

MODEL SCOPE;

     Simulation is of infiltration from a strip source of finite width
located on the soil surface.  Analysis is two-dimensional.  Parameters
included in simulation are time distributed soil-water flux and potential.

INPUT REQUIREMENTS;

     The user specifies the increments in a trapezoidal integration formula,
tolerance for convergence, and strip width as cards in the program.  As
input, the program needs the number of dimensionless x, y, z, and T values at
which solutions are to be determined, the number lines in the analysis, and
the values of the dimensionless variables.  Each parameter is an array
defining values of the parameters at grid points.  Each dimensionless param-
eter is determined from the soil moisture characteristics  (hydraulic conduc-
tivity, moisture content, and pressure head relationships).

SPATIAL AND TIME SCALES;

     Simulation interval is specified by the user although it would be normal
to model a single irrigation event.  The two-dimensional nature, of the model
provides a realistic examination of water distributions over an irrigated
field.

COMPUTER CODE STRUCTURE;

     The computer program utilizes a main program and function routines with
both common and call data transfer.  Program is written in Fortran IV.
Central memory requirements are 24 k bytes and execution time is on the order
of 5-10 central processor seconds per analysis.  Program is readily adaptable
to other computer systems.

BASIC MATHEMATICAL APPROACH;

     The generally accepted moisture flow equation for the time dependent
case becomes a linear differential equation when it is assumed that the
change in moisture content with matric flux potential is linearly related to
the changes in hydraulic conductivity and moisture content.  This assumption
is realistic only when moisture content variations are small such as under
high frequency irrigation systems.  Applications to the soil system are
assumed to be at a steady rate over the interval of simulation.

     The solution to the linearized, time-dependent moisture flow equations
are analytical and presented in the cited reference.   (See model A-4.)
                                      47

-------
GENERAL:
     Computer code must be obtained from the authors.
                                    48

-------
MODEL;                                A-9

REFERENCES;

Wind, G. P., and W. Van Doorne.  1975.  A numerical model for the simulation
of unsaturated vertical flow of moisture in soils.  J. of Hydrology, 24:1-20.

MODEL SCOPE:
     This model simulates the non-steady unsaturated vertical flow in soils.
The model utilizes flow velocity rather than flux to simulate the time and
spatial distribution of water in the profile down to a water table.  Outputs
include the infiltration, vertical movement, water table depths and flux.
Soil conditions at each depth increment are also printed.  Root extraction
and evaporation are not included.

INPUT REQUIREMENTS;

     The user must specify the initial soil moisture, groundwater and surface
conditions, soil hydraulic properties  (conductivity versus moisture tension
and moisture content versus moisture tension), irrigation or rainfall charac-
teristics, geometry of the system  (layers and time increments), and duration
of the simulation.  A number of problems can be evaluated during the same
run.

TIME AND SPATIAL SCALES:
     This model simulates a one-dimensional view of the unsaturated zone.
Time increments of fractions of a day are  integrated over the irrigation or
rainfall interval.

COMPUTER CODE STRUCTURE:
     Program is written  in FORTRAN for a CDC 6600 computer.  Formulas for
storage and execution  times  are given in the reference.  Code consists of a
main program and  subroutine.  Data is transferred with common storage blocks.
Input is accomplished  with cards.

BASIC MATHEMATICAL APPROACH:
     The basic mathematical  expressions utilized by the model are:
 (1) Darcy's  equation  for vertical  flow velocity;  (2) exponential relationship
between hydraulic  conductivity and moisture tension;  (3) conservation of
mass; and  (4) Hooghoudt drainage formula.  The computation procedure begins
with the existing  soil  moisture distribution and the surface input conditions.
The soil is  divided into compartments and  the moisture contents and hydraulic
conductivities for each defined.   Based on averages between compartments  (to
determine  gradients)  the intercompartmental flux is computed.  The compart-
mental moisture  content is then redefined  based on net input-output and the
procedure  is repeated for the next time increment.  Deep percolation is the
flow passing the last unsaturated  boundary.  Infiltration is the flow passing
the first  boundary.

                                      49

-------
GENERAL;

     The computer program is given in the reference cited above.  No
verifications are actually reported, although sensitivity and accuracy are
discussed.
                                       50

-------
                                 APPENDIX B

              DESCRIPTION OF UNSATURATED SOIL SOLUTE CHEMISTRY
                            AND TRANSPORT MODELS
MODEL;                                B-l

REFERENCES;

Dugard, J. D., and M. Reeves.  1976.  Material transport through porous
media:  a finite-element Galerkin model.  Environmental Sciences Division
Publication 733, Oak Ridge National Laboratory, Oak Ridge, Tennessee, March.
201 pp.

MODEL SCOPE;

     Simulation is of transient flow of a dissolved constituent through
either saturated or unsaturated porous media.  Transport mechanisms include
advective transport and hydrodynamic dispersion.  Reactions with the porous
media include chemical adsorption and radioactive decay (if applicable).  The
output describes the time and spatial distribution of the solutes.  The
moisture flow simulation must be handled with another model.  The authors
refer to such a model but it is not incorporated in this work.

INPUT REQUIREMENTS;

     Data must be supplied to define the geometry and initial conditions of
the problem.  Also the properties of the media and solute require definition.
The soil moisture flux is supplied from the output of another program with
tape storage.  Some of the important soil properties are tortuosity, bulk
density, porosity and dispersion coefficient.  The characteristics of the
material being transported include its molecular diffusivity parameters,
radioactive decay rates as applicable and a distribution coefficient for
the linear adsorption function.

TIME AND SPATIALJ3CALES;

     The problem simulated by this model is two-dimensional and of the size
encompased by a field-stream system.  Time increments are set in the program
for calculation purposes, but the period being simulated can extend from
hours to many years.
                                      51

-------
COMPUTER CODE STRUCTURE;

     The program consists of 15 subprograms written in Fortran IV and using
common and call data transfer.  The code is comparatively large (about 2,000
lines) but does not appear to use more than 100-150 k bytes.  Execution times
are dependent on the interval being simulated.

BASIC MATHEMATICAL APPROACH;

     With the water flow phase being given as input, the equations of flow
describe the mass flux of the solute.  Flux is divided into diffusive flux
and advective flux.  Diffusive flux is calculated with a hydrodynamic disper-
sion tensor and advective flux is equal to the Darcy flux of the water times
its concentration.  Adsorption is based on a linear function by assuming the
ratio of per unit mass of solid adsorption to the concentration is constant.
A retardation factor is defined as a linear relationship to the adsorption
constant to delay breakthrough of the constituent.  Some simplification in
the mathematics is achieved by assuming the soil is isotropic with respect to
dispersivity.

     The integration of the differential equations in the spatial frame is
accomplished with the Galerkin finite-element technique using a weighted-
residual method.  Some of the boundary conditions are of the seepage face
and Dirichlet type.

GENERAL;

     The program is demonstrated for the case of seepage return flow to a
stream.  No verification is reported although the numerical solution agreed
well with an analytic solution to the problem studied.  The program is
contained in the reference.
                                     52

-------
MODEL;                                 B-2

REFERENCES:

Jurinak, J.  J., S. Lai, and J. J. Hassett.  1973.  Cation transport in soils
and factors affecting soil carbonate solubility.  EPA-R2-73-235, Office of
Research and Monitoring, U.S. Environmental Protection Agency, Washington,
D.C.  May.  87 pp.

MODEL SCOPE:

     Simulation is of cation transport in soils including processes of
miscible displacement and cation exchange (both linear and nonlinear exchange
isotherms).   Output includes the time and depth distribution of solute and
exchange cation concentrations and an input dump.  Plow is assumed steady.
The nonlinear exchange function is a "Kielland" type.

INPUT REQUIREMENTS;

     User specifies an alphanumeric identification of the analyses, depth and
time increments, total number of time and depth increments, and output
frequency control parameters.  If the Kielland function is used, the input
data set should also include two coefficients for the function.  If the
linear isotherm is used, two other coefficients are needed (slope and inter-
cept) .   The model must be supplied with the dispersion coefficient, median
flow velocity, bulk density, pore fraction, total cation concentration in
solution and the exchange capacity.

TIME AND SPATIAL SCALES;

     This model is a one-dimensional analysis of cation transport in soils.
The time scale for calculations would be on the order of minutes with output
in hours.

COMPUTER CODE STRUCTURE:
     Two similar codes are presented in the reference cited above (one for
each exchange function).  Written for the UNIVAL 1108 computer in Fortran IV,
the codes are very small  (probably less than 10 k bytes) with execution times
ranging from a few seconds to several minutes depending on the time increment
control parameters.  Programs can be easily adapted.

BASIC MATHEMATICAL APPROACH:
     A material balance equation is written in an explicit finite difference
form to describe the exchange and miscible displacement of a steady downward
flux.  Numerical dispersion is eliminated with proper grid spacing.

GENERAL:
     The computer programs are contained in the reference.  Verification is
reported with laboratory data.

                                      53

-------
MODEL;                                B-3

REFERENCES;

Oster, J. D., and B. L. McNeal.  1971.  Computation of soil solution
composition variation with water content for desaturated soils.  Soil Sci.
Soc. Amer. Proc., 35:436-442.

Oster, J. D., and J. D. Rhoades.  1975.  Calculated drainage water composition
and salt burden resulting from irrigation with river waters in the western
United States.  J. Environ. Qua!., 4:73-79.

MODEL SCOPE:
     This simulation involves the salinity composition of root zone drainage
flows.  Parameters include dissociated ionic species, gypsum and lime dissolu-
tion/precipitation, SAR, CO- partial pressure at bottom of root zone  (Pco ),
pH, EC, TDS, ionic strength, activities, and activity coefficients.  Output
also includes CO. and HCO~ ion pairing, and 1st and 2nd dissociation constants
for carbonic acid.  The computer code features an input dump, error flags for
iterative sections, but no intermediate results.
INPUT REQUIREMENTS;

     Control parameters are provided so a user may include magnesite  (not
fully developed) and forced saturation of soil lime.  Up to 100 individual
analyses can be run without dimension and common modification.  Data include
an identification of the test, Pc©2 at lower root zone boundary, pH, Na+,
Mg++, Ca++, K+, HCO-j, Cl~, and SO^ in applied irrigation water, and the
estimated leaching fraction.

SPATIAL AND TIME SCALES:
     Simulation is of a single or representative site in an irrigated field.
Analysis is steady-state and encompasses the drainage period of the irriga-
tion interval.  Model is one dimensional.

COMPUTER CODE STRUCTURE;
     Internal program linkage is accomplished by subroutine and functions
with common data transfer.  The program can be easily adapted, all or in
part, to other computer systems or programs.  Approximate core and execution
time requirements are 70 k bytes and 10-15 central processor seconds,
respectively.  Programmed originally for an IBM 360 computer using a Fortran
IV language.

BASIC MATHEMATICAL APPROACH:
     The basic modeling approach along with important assumptions, supportive
literature, and data requirements is well described in the previously cited
references.  The interested user is therefore referred directly to these
sources for more detailed information.  In general terms, the computations

                                      54

-------
begin by pre-concentrating the ions in the irrigation water by dividing each
by the leaching fraction.  Then an iterative procedure  (governed by con-
vergence tolerances on ionic charge balance, changes in ionic concentrations
for Ca and CO^, and dissociation constants) is used to compute the salinity
constituents in the drainage.  The extended form of the Debye-Huckel equation
is utilized to calculate activity coefficients in the equilibria analyses.
Gypsum and lime precipitation/dissolution is used to achieve convergence.
Leachate pH is governed by the soil lime equilibria and cation exchange is
not included.

GENERAL:
     The computer code is available from the authors cited above.
Documentation within the code is excellent and substantial verification
trials are reported to demonstrate the utility and limitations of the model.
                                     55

-------
MODEL;                                B-4

REFERENCES;

Rai, D., and W.T. Franklin, 1973.  Program for computing equilibrium solution
composition in CaCCXj and CaSC>4 systems from irrigation water compositions.
Water Management Technical Report No. 29.  Colorado State University,
Fort Collins, Colorado.  October.  42 pp.

MODEL SCOPE:

     Simulation is of the salinity composition of root zone drainage.
Parameters include dissociated ionic species, gypsum and lime dissolution -
precipitation, SAR, pH, EC, TDS, ionic strength,_activity, and activity
coefficients.  Output includes CO^, HCO3, and 804 ion pairing with Ca   ,
Mg++, K+, and Na"*~.  The computer code features an input dump, error flags for
iterative sections, but no intermediate results.

INPUT REQUIREMENTS;

     The program does not include any control options.  Input data are sample
identification, total concentrations of Ca"1"1", Mg++, Na+, K+, HCO^, 504, and
Cl~ in  the irrigation water, CO_ partial pressure, and the leaching fraction.
Units are in equivalents per liter except for CO2  (atm) and leaching fraction
 (expressed as a fraction).  Program operates on one sample during each run.

SPATIAL AND TIME SCALES;

     The simulation involves a single or average condition in an irrigated
field on a one-dimensional basis and steady-state system is assumed.

COMPUTER CODE STRUCTURE;

     The code consists of only a main program in the latest version, although
the program in reference cited above contained a subroutine for solution of
cubic roots.  The program can be easily modified for use on computers other
than the CDC 6400 at Colorado State University.  Core requirements are
approximately 30 k bytes.  Execution time per run is 15-70 central processor
seconds.  Code language is Fortran IV.

BASIC MATHEMATICAL APPROACH;

     The basic modeling approach, along with pertinent assumptions and
simplifications, is well documented in Rai and Franklin (1973).  Program is
now revised to include nesquehonite.  Computations begin by preconcentrating
ionic concentrations in the irrigation water by dividing each by the leaching
fraction.  Ionic strength and activity coefficients are initially calculated
by assuming no pairing.  Then, a three-dimensional successive approximation
loop is entered in which calcium equilibrium is used to arrive at leaching
water chemistry.  The primary loop flow involves (1) computation of free
anions;  (2) calculation of free cations; (3) calculation of charged ion
pairs;  (4)  calculation of ionic strength; and (5) calculation of activity

                                      56

-------
coefficients.   Then,  lime or gypsum equilibrium is checked.  If solution is
in equilibrium, the program goes to next calcium system or concludes.  If
nonequilibrium exists, approximations resume at beginning of loop.  The
Debye-Huckle equation is used for activity coefficients and literature
references for equilibria constants are used for the calcium system.  Cation
exchange is not included.

GENERAL;

     The code is available from the authors cited above.  Code is moderately
documented itself whereas the reference is extensive.  Model has been verified
against literature and field data.  This program is also available for an
HP 9825 programmable calculator from the writer.
                                      57

-------
MODEL:                                B-5

REFERENCES:

Saxton, K.E., G.E. Schuman, and R.E. Burwell.  1977.  Modeling nitrate
movement and dissipation in fertilized soils.  Soil Sci. Soc. Amer. Proc.,
41:265-271.

MODEL SCOPE;

     Simulation is of the nitrate-nitrogen-occurrence, movement, and
dissipation in 1-2 meters of the unsaturated zone under fertilized agricul-
tural lands.  Processes included in the model are solute transport  (infil-
tration, redistribution, and percolation), crop uptake, disposition, of added
fertilizer, rainfall additions, mineralization and nitrification.  Denitri-
fication, fixation, and runoff losses were assumed negligible.  Other nitrogen
processes were considered negligible.  The moisture flow requirements must be
input to the model and must include soil evaporation, plant adsorption from
each 15 cm layer, infiltration, soil moisture redistribution between layers,
and soil moisture in each layer.

INPUT REQUIREMENTS;

     Water flow data for this model comes from the model described by Saxton,
et al.  (1974).  Output from that model or similar daily values for uptake by
soil segment, infiltration, soil moisture storage volumes, and deep percola-
tion must be incorporated.  Then, the time distributed data describing
nitrogen additions to the system, total and depth distributed plant uptake
of nitrogen, and the initial nitrogen profile in the soil are read.

TIME AND SPATIAL SCALES:
     Daily calculations are made of a one year event in a one-dimensional,
transient, unsaturated zone.

COMPUTER CODE STRUCTURE:
     Using a main program to read tape output from the soil moisture model,
the general nutrient calculations are carried out in subroutines using common
and call statement data transfer.  Execution times for the Fortran Code are
on the order of 10-20 seconds.  Core requirements should not exceed 50-60 k
bytes.

BASIC MATHEMATICAL APPROACH:
     Nitrogen transport is based on the assumption that within each 15 cm
soil segment, the concentrations of nitrate are uniformly distributed.  Thus,
nitrates leaving or entering a layer are calculated as the product of the
moisture movement and the existing nitrate concentrations.  New concentra-
tions are computed for each hour and soil layer to aggregate the transient
system on a daily basis.  The withdrawal of nitrates from the soil layers is
assumed proportional to water absorption.  An annual uptake figure is

                                     58

-------
distributed on a daily basis assuming a proportionality with dry matter
production and then distributed in the soil profile according to root
distributions.  Nitrogen added to the upper soil layers in precipitation or
fertilizer moves as indicated above, but ammonium fertilizers are assumed to
be either adsorbed or converted to nitrate.  A  first order temperature
dependent decay function was used to describe the ammonium conversion to
nitrate.  A 10% loss due to volatization was assumed on initial application.
Mineralization is fixed at 52 kg/ha per year within the upper soil layer and
is distributed with time proportionally to the temperature in the soil
between May and September.

     The calculation proceeded in the following sequence:  fertilizer
addition, mineralization addition, uptake, infiltration addition and
transport, and transport by redistribution and percolation.  A revised
nitrate profile is calculated after each process and printed at the day's
end.

GENERAL;

     The model has been verified  (with good agreement) with field data from
two watersheds.  This is a readily usable model only if the moisture flow is
supplied by the earlier model.  Thus, modifications should be made to operate
this one independently for widespread application.  Computer codes must be
requested from the authors of the cited references.
                                      59

-------
MODEL^                                B-6

REFERENCES:

Tanji, K. K., L. D. Doneen, G. V. Ferry, and R. S. Ayars.  1972. Computer
simulation analysis on reclamation of salt-affected soils in San Joaquin
Valley, California.  Soil Sci. Soc. Amer. Proc. 36:127-133.

MODEL SCOPE;
                                                           I _t_    _i	i     i    ^
  _  Simulation is of the interaction of soluble salts  (Ca  , Mg  , Na , K ,
304, Cl~, HC03, and N03) and Boron in a soil profile.  The model was developed
to describe leaching processes so the output includes drainage water quality,
soil profile chemistry  (changes in soluble, and adsorbed salts).  It also
predicts the depth of leaching water required to achieve specific values of
soluble salts, boron, and exchangeable sodium concentrations in the soil
profile.  Output includes ionic strengths, activity coefficients, ion pairs,
stoichiometric solute concentrations  (dissociated and undissociated), total
soluble cations, SAR, and ESP.

INPUT REQUIREMENTS;

     Input data include leaching water salinity composition, initieil soil
profile salt distribution, depth distributed exchangeable Ca  , Mg  , Na .
cation exchange constants, gypsum content, adsorbed boron, Langmeier
constants, and moisture contents at saturation and field capacity.

TIME AND SPATIAL SCALES;

     The principal time frame is the period covered by an application of
irrigation water and the subsequent leaching.  The model is one-dimensional
and encompasses the unsaturated zone of depth specified by the user.

COMPUTER CODE STRUCTURE;

     The program is a Fortran code without subroutines or user supplied
functions.  Core and execution times are not given but appear small  (40-60 k
bytes, and about 1-2 minutes).

BASIC MATHEMATICAL APPROACH;

     The soil profile is divided into layers with moisture movement between
layers when field capacity is exceeded.  A chromatographic approach to salt
transport between layers is assumed.  Solute from one layer moves into a
lower layer, mixes with the resident solution (at equilibrium) and
equilibrates.

     The calculation of the equilibrium chemistry for the salts utilize the
Guggenheim-Davis equation for activity coefficients and experimental
solubility and dissociation coefficients.  Boron adsorption-desorption is
simulated with a Langmeier isotherm.
                                      60

-------
GENERAL:
     Model was verified with field data and has served as a starting point
for several steady-state models developed to predict the chemistry of
leaching water.  Code must be requested from authors cited above.
                                     61

-------
                                APPENDIX C

              DESCRIPTION OF UNSATURATED ZONE HYDROLOGY MODELS
MODEL:                                C-l
REFERENCES:
Feddes, R. A. and H. Zaraday.  1977.  Numerical model for transient water
flow in nonhomogeneous soil-root systems with groundwater influence.
Proceedings of IFIP Working Conference on Modeling and Simulation of
Land, Air, and Water Resource Systems, Ghent.  September.  18 pp.

MODEL SCOPE:
     Simulation is of the unsaturated nonhomogeneous soil system under
transient input conditions and root extraction.  The uptake of water by the
root system  is incorporated  into the one-dimensional moisture flow equation
as a sink term dependent on  soil moisture head, effective rooting depth, and
potential transpiration.  Output includes the time and depth distributed soil
moisture conditions, the evapotranspiration, and the interactions with a
water table.

INPUT REQUIREMENTS;

     The soil is divided into 20 increments and described with the following
data:  pressure head—hydraulic conductivity—soil moisture content relations;
root depth with time; a time distributed root depth correction to account for
nonactive roots; pressure head at the wilting point, anaerobiosis point, and
at the optimal uptake level; initial soil moisture and water table conditions;
and potential transpiration.

TIME AND SPATIAL SCALES:
     This model is a one-dimensional analysis of the unsaturated zone.
Calculations simulate events over a daily period of as long as desired (one
year in reference).

COMPUTER CODE STRUCTURE:
     The code is a Fortran IV program consisting of a main program and a
plotting subroutine.  Neither core nor execution times are.Riven but should
not be excessive for most computer systems.
                                      62

-------
BASIC MATHEMATICAL APPROACH;

     The one-dimensional soil moisture flow equation is solved by an implicit
finite difference technique.  The root extraction sink term is essentially
the same as described by Neuman et al. (1975) based on earlier work by
R. A. Feddes, although slight modifications are made to the effective root
distribution.  The effect of the water table is to provide upward moisture
movement in response to root extractions and downward movement due to deep
percolation.

GENERAL:
     This is a simplified version of the model described by Neuman, et al.
 (1975) and Ostrowski  (1976) but verification against field data showed better
accuracy.  Code must be requested from authors cited in reference.
                                      63

-------
MODEL;                                C-2

REFERENCES:

Gupta, S.K., K.K. Tanji, D.R. Nielsen, J.W. Biggar, S.C. Simmons, and
J.L. Maclntyre.  1977.  Field simulation of soil-water movement with water
extractions.  Water Science and Engineering Paper 4013.  University of
California, Davis, California.  May.  95 pp.

MODEL SCOPE;

     Simulation of "infiltration, vertical seepage, and uptake by plants as
related to the hydraulic properties of the soil, soil layering, the root
growth characteristics of a given crop in a given soil profile, evapotrans-
piration rates, and frequency, rate, and amount of irrigation," is presented.
Outputs, in addition to an input dump, includes hourly values of moisture
content, suction, and flux throughout the profile during an irrigation and
daily values during nonirrigation periods.  Daily mass balances are also
given along with distributed totals for each run describing the irrigation
input, evaporation and transpiration, and flux at each node.  Data may be
supplied by disk or cards.

INPUT REQUIREMENTS;

     The program requires a number of control parameters to direct the
calculations through the alternatives for various steps.  First, boundary
conditions when not calculated directly must be specified for field and crop
conditions.  The problem size and decomposition structure is defined next,
i.e., distribution of nodes and material numbers above and below each node.
Dates for planting and harvesting are also read in; then a series of flags
are specified.  Leaf-area-index options include direct input of LAI versus
time, or a user specified distribution.  Partitioning transpiration and
evaporation from inputed values of potential evapotranspiration can be either
computed with techniques using LAI related energy balance at the soil surface,
or an option supplied by the user.  Stress effects on transpiration include
options for on-off, logarithmic decrease, linear decrease, a combination on-
off and linear or a user supplied function.  Soil profiles may be either
homogeneous, nonhomogeneous, or layered, each of which requires the hydraulic
conductivity-moisture content relationship (expressed as polynomials).  Root
growth options involving a negative exponential relationship, distributed
density functions, or a user supplied system is specified and associated data
read into the model.  Moisture extracted by the rooting system is evaluated
in a sink term added to the moisture flow equation.  The form of the sink may
be a macroscopic Nimah-Hanks type, related to soil suction, or supplied by
the user.  Associated data would then be necessary to define the pertinent
coefficients.  Finally, surface boundary conditions describing static,
quasi-dynamic, dynamic, or semi-infinite depth must be known.
                                      64

-------
TIME AND SPATIAL SCALES;

     During an irrigation, time increments may be as small as 0.1 to 1.0
hour with aggregated values presented daily.  Spatial resolution is a one-
dimensional analysis of the unsaturated region below an irrigated surface.

COMPUTER CODE STRUCTURE;

     A main program-subroutine system with common and call information
transfer constitutes the essential code structure.  Code is written in
Fortran IV.  Execution time and core requirements are not specified in the
reference, but neither are expected to be a problem at most computer
facilities.

BASIC MATHEMATICAL APPROACH;

     Water entering the soil via irrigation or precipitation either
redistributes in the root zone or evaporates from the soil surface.  Knowing
the potential evapotranspiration rate, surface evaporation is delineated by a
leaf-area-index related radiant energy balance computation at the soil
surface.  The moisture not evaporating from the soil surface moves through
the soil profile and/or is extracted by the rooting system.

     Soil moisture movement in the root zone is simulated with the single
phase, time dependent solution to the Darcy equation.  A sink term is added
to the basic finite difference solutions to accomplish plant uptake.  The
model handles the sink term in an interesting way.    An iterative solution
of the K-9 system is undertaken until a mass balance between internodal flux
and total transpiration is accomplished.  A similar mass balance relating the
moisture remaining after plant uptake with the distribution of moisture in
the profile is utilized to iteratively change time resolutions until both
estimates are congruent.  Thus, the model automatically insures a mass
balance by forcing the detached solutions to the flow equations to be in
agreement.

GENERAL;

     This model is the first phase of a comprehensive study aimed at
simulating nitrogen behavior in soils.  Verification has been made against
very detailed field data, although only at one location.  The code reflects a
fair state-of-the-art description of the physical processes in the root zone
even though more refined treatments of some segments are possible.  However,
the sophistication among these segments appears quite homogeneous.  Although
classified primarily as a research tool, the approach allows application to a
number of field situations.  Code is available from the authors cited above.
                                     65

-------
MODEL;                                C-3

REFERENCES;

Hillel, D.  1977.  Computer Simulation of Soil-Water Dynamics.  International
Development Research Centre, Ottawa, Canada,  pp. 61-78.

MODEL SCOPE:
     Model simulates non-isothermal evaporation from a soil surface based on
an energy balance approach.  Moisture transport is simulated as a two phase
flow problem and soil heat flux is incorporated.  Output includes a time and
spatial distribution of soil moisture and heat, surface energy balance, and
evaporation.

INPUT REQUIREMENTS:
     The code requires a description of initial soil moisture condition, soil
properties, and atmospheric inputs.  Soil hydraulic properties include the
moisture conductivity-suction relationship, soil thermal conductivity  (air,
water, and solids), volumetric heat capacity  (air, water, and solids),
albedo, field capacity, porosity, surface roughness, soil temperature-water
vapor thermal conductivity function, and the geometry of the system being
modeled.  Atmospheric variables listed as input are time distributed air
temperature, dew point temperature, wind speed, and solar radiation.

TIME AND SPATIAL SCALES:
     The basic unit of time is one hour with periods of several days being
evaluated.  The spatial scale is a one-dimensional evaluation of the first
1-2 meters of the soil profile.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH:

     The surface energy balance includes approximately the same approach
described by the Penman equation for evapotranspiration.  In this case,
the model uses more of the analytical functions and less regression type
relationships.  Water flow occurs according to Darcy's Law and conservation
of mass.  Heat flow is described by Fourier's Law written in one-dimension.

     The modeling approach involves three steps.  First, the hydraulic and
thermal gradients and associated conductivities between layers are determined
based on initial or existing conditions.  Then, the heat, vapor, and water
flux is determined by multiplying the gradients and conductivities by the
time increment.  And finally, the contents of each depth increment or com-
partment are updated by evaluating the net inflows.  The volumetric condition
of each soil segment is determined based on average gradients existing


                                     66

-------
between compartments.  The iterative solution involves a fourth order Runge-
Kutta procedure.

     At the soil surface, both water and heat fluxes are calculated.  Surface
temperature is used to segregate the latent heat flux  (evaporation) and the
sensible heat  (soil warming or cooling of the air).  The net energy for
evaporation at the soil surface and the net vapor transport through the soil
profile are successively evaluated with surface temperature as the dependent
variable until convergence is satisfied.  Then the surface evaporation and
soil heat flux are known for this time interval and the calculation proceeds
through the soil profile.

GENERAL:
     This model is apparently a second or third generation version developed
in recent years by various researchers in conjunction with Dr. Hillel.  In
this reference, the model is applied to a Gilat fine sandy loam soil in
Israel and used to evaluate the effect of albedo variance.  These results are
simulation  runs, not calibration or verification trials.  Code is included
in reference.
                                      67

-------
MODEL;                                C-4

REFERENCES:

Hillel, D.  1977.  Computer Simulation of Soil-Water Dynamics, International
Development Research Centre, Ottawa, Canada,  pp. 79-130.

MODEL SCOPE:
     Simulation is of infiltration, redistribution, evaporation, and deep
percolation in the unsaturated root zone.  Output includes the time and
spatial distribution of water in the system.  Results can be plotted and
input data can also be listed.

INPUT REQUIREMENTS;

     The vertical soil profile is divided into a number of compartments or
layers of variable thickness.  Each layer is assumed to have an initial
moisture content, but variable values can be assigned.  The various layers
can also have different hydraulic conductivities if desired.  Other input
consists of the moisture content-hydraulic conductivity-suction relation-
ships, and the evaporation-rain or irrigation input frequencies.

SPATIAL AND TIME SCALES;
     The spatial scale is the one-dimensional compartmentalization of the
unsaturated zone under consideration.  The governing equations are solved in
intervals of one second with periods of several hours or days considered.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The basic equations of flow are divided into parts.  The flux within the
soil profile is calculated by the CSMP system of iteratively calculating
moisture content, determining average conductivities and gradients, calcu-
lating interlayer flux and recomputing moisture contents based on net inflow-
outflow.  The intensity of rainfall or irrigation is compared to intake rate
of the soil surface to delineate infiltration and surface storage.  Evapora-
tion at the soil surface occurs at the potential rate (input data) until the
upward flux or input is less than this rate.  And finally, the deep percola-
tion is taken as the lowest compartment's hydraulic conductivity.  Uptake and
transpiration are not included.  A mass balance is maintained as a check on
the computations and output summary.

GENERAL:
     Analysis of the model does not indicate either verification or
calibration trials.  Rather, evaluation of soil types and initial
conditions are presented.  Code is included in the reference.

                                      68

-------
MODEL:                                C-5
REFERENCES:
Saxton, K. E., H. P. Johnson, and R. H. Shaw.  1974.  Modeling
evapotranspiration and soil moisture.  Trans. ASAE, 17:673-677.

MODEL SCOPE:
     The simulation is of the evapotranspiration and soil moisture profiles
in an irrigated system.  Evapotranspiration is determined as the sum of
canopy evaporation, soil evaporation, and plant transpiration.  Soil moisture
is redistributed by hydraulic conductivity-tension functions.

INPUT REQUIREMENTS;

     The first input data are the time and depth distributed root distribution
percentages to proportion the root water extraction patterns.  Then, daily
values of potential evapotranspiration are read in based on a suitable
approach (combination equation and pan evaporation in this work) as are net
radiation values.  Then, the data describing the stress curves are read in to
adjust the crop transpiration.  For each soil type present, the user must
detail the conductivity-tension relationship, wilting point, field capacity,
saturation and initial soil moisture.  Water inputs to the system through
infiltration are input along with canopy development curve data.  Output
can be directed to paper, tape, or both types of printout.

SPATIAL AND TIME SCALES:
     The time interval for calculation is one or two irrigation seasons
evaluated a day at a time.  Spatial scale is a one-dimensional profile
through the 1 to 1.5 meter root zone.

COMPUTER CODE STRUCTURE:
     The program is contained primarily in a main program of sequential
calculations.  A series of subroutines using call statement data transfer
exclusively perform interpolation and data adjustment functions.  Code is
written in Fortran IV.  Core requirements noted in printout appear to exceed
150 k bytes, but should not exceed 40-50 k in most computers.  Execution
times are on the order of 10-20 central processor seconds.

BASIC MATHEMATICAL APPROACH:
     The basic model calculations involve three sequential steps performed
repeatedly on a daily interval:  (1) computation of actual evapotranspira-
tion; (2) adding infiltration to soil moisture; and (3) redistributing the
soil moisture.  Beginning first with input values of potential evapotranspi-
ration,  the initial calculations segregate canopy evaporation from soil
evaporation and plant transpiration.  Maximum interception canopy evaporation
is set at 2.5 mm.  Potential values for soil evaporation  (assumed to come
only from first 15 cm) and transpiration are divided according to the canopy

                                     69

-------
shading percentage.  Reduction in potential soil evaporation to an actual
value is made with a literature ratio expressed as a function of soil
moisture.  Transpiration is adjusted by crop curves.  The daily transpiration
(and evaporation in the first 15 cm) is assigned to respective soil segments
based on root distributions.  These values are then adjusted further for
stress using literature relationships.  After the root extraction has been
deducted from the soil moisture, infiltration due to daily precipitation (or
irrigation if present) is evaluated.  All infiltration up to 90% of the
available capacity in each soil segment is initially stored therein beginning
with the top layer.  And finally, a four hour steady-state analysis of the
one-dimensional Darcy equation for unsaturated flow utilizing the conductivity-
tension relationship is used to redistribute the moisture in the profile for
the next day's analysis.

GENERAL:
     Three years of detailed field data were used to verify the modeling
approach.  Calibrations showed good agreement with measured data.  This model
does not include extensive treatment of each individual process which would
make it exclusively a research model.  The code is available from the authors
cited.
                                     70

-------
MODEL;                                C-6

REFERENCES;

Trava, T. L.  1975.  Allocating available irrigation water on the farm.
Unpublished Ph.D. Thesis, Department of Agricultural Engineering, Colorado
State University, Fort Collins, Colorado.  December.  107 pp.

Jensen, M.E., C.E. Franzoy, and D.C.N. Robb.  1970.  Scheduling irrigations
using climate-crop-soil data.  J. Irr. and Drain. Div., ASCE, 96:25-38.

Kincaid, D.C. and D.F. Heerman, 1974.  Scheduling irrigations using a
programmable calculator.  ARS-NC-12.  Agricultural Research Service, U.S.
Department of Agriculture.  February.  55 pp.

MODEL SCOPE;

     Simulation of irrigated crop evapotranspiration is based on climatic
data utilized in a combination energy flux-advection model (Modified Penman
Equation).  Parameters in the analysis include the potential evapotranspira-
tion for 10 field crops, actual use based on stress and growth stage adjust-
ments, predicted dates of future irrigation, and soil moisture status
estimates.

INPUT REQUIREMENTS;

     The program is set up to compute multiple crop, multiple farm, multiple
region irrigation updates and schedules.  After this problem structure is
defined, the user supplies a region by region,function for the average
expected evapotranspiration.  Then soil moisture, consumptive use, climatic,
and irrigation data for the three previous days are read in and evaluated.
Crop data are also read in with each run.

SPATIAL AND TIME SCALE;

     Evapotranspiration and irrigation scheduling is an areal average,
nonspatial computation.  Time scale is on a daily basis, and the model is
one-dimensional.

COMPUTER CODE STRUCTURE;

     Main program and subroutine structure, written in Fortran IV, compromise
the Trava version of the model.  Common data transfer is primarily employed.
Core requirements amount to about 100 k bytes but some of this involves the
optimization package incorporated in the Trava model.  Execution time is on
the order of 50 central processor seconds.

BASIC MATHEMATICAL APPROACH;

     During the course of a scheduling run, detailed soil, crop, and climatic
data for the period since the last update are used to update the present time
frame.  The modified Penman Equation is used as the basic tool with stress,
recent irrigation, and crop corrections used as described by Jensen, et al.

                                      71

-------
(1970) and Kincaid and Heerman (1974).  The future prediction utilizes the
same approach with the exception that a mean expected evapotranspiration
function is used instead of the Penman equation.

GENERAL;

     The Penman equation is used extensively in the more detailed simulation
models and this model's approach to calculating consumptive use is among the
most correct.  Thus, although it is more of an operational tool, it is
included in this writing because of its widespread applicability and basis
for consumptive use calculations in other models.  Code is only available
from the writer.
                                      72

-------
MODEL:                                C-7
REFERENCES:
Neuman, S. P., R. A. Feddes, and E. Bresler.  1975.  Finite element analysis
of two-dimensional flow in soils considering water uptake by roots:  I.
Theory, II. Field Applications, Soil Sci. Soc. Amer. Proc., 39:224-237.

Ostrowski, J. T.  1976.  Behavior of groundwater subject to irrigation of
effluent—a case study.  M.S. Thesis, Civil Engineering Department, University
of Maryland, College Park, Maryland.  121 pp.

MODEL SCOPE:
     Simulation is of infiltration, evapotranspiration, root extraction, and
soil moisture distribution in an irrigated field.  Soil profiles extend to a
water table lower boundary.

INPUT REQUIREMENTS:

     The geometry of the soil system is first defined by locating a finite
element grid.  Then the soil hydraulic properties are incorporated (porosities,
permeabilities, tension versus moisture content relations, and relative
permeability versus moisture content relations).  Initial and boundary
conditions at each element node require description of potential evaporation
and transpiration, infiltration, and flux.  And finally, the time limits of
each problem are set.

TIME AND SPATIAL SCALES:
     The model is a two-dimensional analysis of an irrigated field condition.
Calculations are made for a period of hours or parts of hours and extend as
long as desired.

COMPUTER CODE STRUCTURE:
     The code as modified by Ostrowski  (1976) is written in Fortran IV.  A
main program with subroutines  (common and call data transfer) form the basic
structure.  Execution times are very large, one week, for example, requires
20 mintues.  Core requirements are not noted because of tape storage for
data, however, the program should be useable on computers with available core
in excess of 50 k bytes.

BASIC MATHEMATICAL APPROACH:
     The basic soil moisture flow equations are solved using the Galerkin
finite element method, using a discretization scheme.  Two aspects of the
flow equations merit attention here.  Potential evapotranspiration in the
Ostrowski (1976) version is based on pan evaporation adjusted by a tempera-
ture dependent growth index and the porosity-tension characteristics of the
soil.  Surface evaporation is a prescribed fraction of the total evapotran-
spiration with potential transpiration being the difference.  The uptake of


                                      73

-------
moisture in a soil segment is assumed to be a function of hydraulic
conductivity multiplied by the difference in pressure head difference between
the root-soil interface and the surrounding soil and divided by a propor-
tionality coefficient describing the exponentially distributed root mass.
Actual transpiration is calculated as the maximum sum of water extractions in
each soil segment subject to a limit imposed by the potential rate.  Evapora-
tion is only allowed to occur in the first soil segment.  The flux of water
at the soil surface (after evaporation) is an input variable.

GENERAL:
     The model in its various forms has been verified against field data at
several locations.  The practical use of the program, however, is limited by
available computer time and the complexity of required data.  The code is
available from the authors cited above.
                                     74

-------
                                APPENDIX D

            DESCRIPTION OF UNSATURATED ZONE HYDRO-QUALITY MODELS


MODEL;                                D-l

REFERENCES:
Beek, J. and M. j. Frissel.  1973.  Simulation of nitrogen behavior in soils.
Centre for Agricultural Publishing and Documentation, Wageningen, The
Netherlands.  67 pp.

MODEL SCOPE:
     Simulation is of the nitrogen transformation and transport in cultivated
soils under aerobic conditions and having a large buffer capacity.  Transport
and temperature are modeled with the models described by DeWit and Van Keulen
(1975).  The analysis involves the microbial decomposition of organic matter,
and subsequent processes of mineralization, immobilization, humus formation,
nitrification and leaching.  Denitrification and the soil pH reactions are
not included.  Output includes the time and .depth distribution of the nitrogen
compounds and the removal of nitrogen from the system through leaching.  The
soil profile is divided into 20 compartments or layers.  Nitrogen transforma-
tions occur in only the first 4 layers, leaching in all 20, and heat flow in
the first ten.  Initial conditions limit water content, temperature, and
solute characteristics to the same in each layer.

INPUT REQUIREMENTS;

     The input data to this model are as follows:  layer thickness, initial
moisture content, initial nitrate-nitrogen content, initial temperature, heat
capacity (water and minerals), relative amount of minerals, initial mass of
proteins, sugars, lignin, cellulose and biomass, initial ammonium, number of
ammonium oxidizers per soil volume, maximum nitrification rates, bulk density,
hydraulic conductivity and diffusivity as a function of moisture content,
saturation values, water content at lower boundary, irrigation frequency,
irrigation depths, application rates, beginning of irrigation, irrigation
durations, maximum hydraulic head at surface, heat conductivities versus
moisture contents, time distributions of surface temperature, dispersion and
diffusion coefficients, tortuosity, carbon dioxide parameters, and the array
of nitrogen transformation coefficients, and decomposition rates.
                                      75

-------
SPATIAL AND TIME SCALES;

     The basic time increment is the irrigation interval evaluated by
calculation intervals down to minutes.  The simulation is a one-dimensional
treatment of the first 1-2 meter crop root zone area.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     Since the water, heat, and mass transport processes are described in
this work under the DeWit and Van Keulen  (1975) model, this discussion will
be limited to the nitrogen system.

     Organic matter is assumed to consist of protein, sugar, cellulose, and
lignin which resist microbial attack in descending order.  Decomposition
rates are not dependent on concentrations but are both temperature and
moisture dependent.  Humus and biomass also represent organic matter; humus
the nondecomposable form.  The decomposition is a first-order kinetic
approach proportional to the amounts of each segment of organic matter
available.  Temperature and moisture contents other than optimal values are
adjusted by proportionality constants.  The use of ammonium and nitrate by
the biomass is based on carbon/nitrogen ratios and are related by a first-
order function of the decomposition rates.

     The nitrification segment of the model assumes that nitrate moves with
the moisture but ammonium is adsorbed completely  (NH^ exchange processes are
considered).  The rate of NH^ to N03 conversion is determined as the ammonium
to nitrate ratio with other reactions occurring faster and, therefore, are
ignored.  This rate is proportional to the number of ammonium oxidizers and
is independent of the ammonium present.  Both suboptimal temperatures and
moisture content are handled as in the previous paragraph.

GENERAL:
     Verification is not discussed.  A computer code listing is contained in
the reference.
                                     76

-------
MODEL;                                D-2

REFERENCES;

Davidson, J.M., G.A. Brusewitz, D.R. Baker, and A.L. Wood.  1975.  Use of
soil parameters for describing pesticide movement through soil.  EPA-660/2-
75-009.  National Environmental Research Center, Corvallis, Oregon.  May. 162
pp.

MODEL SCOPE;

     This model simulates the simultaneous transport of water and solutes
through a field soil.  The analysis describes the transient and interactive
nature of the solute-soil matrix system.  Assuming isothermal conditions, the
simulation evaluates pesticide concentrations in both the solute and soil
matrix with time and depth.  The transport analysis is a miscible displace-
ment approach including convection, diffusion, mechanical dispersion, and
adsorption-desorption.  Output from the program includes a partial listing of
input control parameters and data, sorbed and solute concentrations as
functions of time and soil depth, and intermediate values of coefficients
and water-solute flux.  Controls are included for iterative calculations.

INPUT REQUIREMENTS;

     In addition to a number of control parameters setting up the alternative
program options, the user must supply a comprehensive list of soil hydraulic
properties, boundary conditions and chemical concentrations and character-
istics.  A series of optional data sets are necessary depending on the
structure of the program selected.

     Soil physical and hydraulic data needed include:  soil depth, saturated
hydraulic conductivity and moisture content, soil evaporation rates and
duration, moisture content versus pressure head relationships, moisture
content versus hydraulic conductivity data, and bulk density.  Boundary
conditions include time and depth increments, time and depth limits, initial
soil moisture conditions, and characteristics of water applications at the
surface.  Other input parameters are included in the description of the
boundary conditions which establish certain system reactions, but the reader
is referred to the reference material for more specific details.  The chemical
data which are needed in the program consist of dispersion, adsorption,
desorption and vapor diffusion coefficients, concentration of pesticide or
chloride in the infiltrating moisture, and initial conditions within the soil
profile.

SPATIAL AND TIME SCALES;

     This simulation is a one-dimensional analysis of pesticide transport
through an unsaturated soil profile encompassing a crop root zone and the
strata below.  Time intervals for intermediate calculations are set by
program user but overall simulation encompasses the application-drainage
periods.


                                      77

-------
COMPUTER CODE STRUCTURE;

     A main program-subroutine format is used with call and common data
transfer.  Code is written in Fortran IV.  Listing is well documented with
comment cards for easy application of the program.  Input and output require
no special tape or disk apparatus.  Model should be easily adapted to computer
systems other than the IBM 360 originally used.

     Core and execution times are not noted in reference material but do not
appear to be limitations for facilities having in excess of 100 k internal
memory and allowing several minutes in central processor execution time.

BASIC MATHEMATICAL APPROACH;

     The analysis begins by using a Crank-Nicolson implicit finite difference
scheme to compute the change in water content as a function of time through-
out the soil profile.  The moisture flow equation, based on single phase flow
and derived from Darcy's equation and mass balance, is evaluated in the soil
matrix by solving its analog form (pressure head) and then calculating
moisture content using head-moisture content relationships.

     The chemical aspects of the problem begin with assumptions that the
adsorption-desorption process can be adequately simulated by single-valued
functions.  Linear relationships are assumed between adsorbed and solute
pesticide concentrations when a pesticide-soil characteristic, N, equals 1
(Freundlich Equation).  The rate of mass transfer to the adsorbed phase is
determined with rate coefficients based on the Arrhenius Equation.

     The transport of pesticides in the soil system under transient conditions
is described by a miscible displacement equation which includes expressions
for convection, diffusion, mechanical dispersion, and adsorption-desorption.
The resulting differential equation is solved with an explicit finite
difference procedure with a number of steps taken to correct any numerical
dispersion that may occur.

GENERAL;

     This model has been tested under both laboratory and field conditions.
The results appear sufficient to suggest use of the model.  The code itself
is contained along with its documentation in the reference cited above.
                                     78

-------
MODEL;                                D-3

REFERENCES:

Davidson, J.M., G.A. Brusewitz, D.R. Baker, and A.L. Wood. 1975.  Use of soil
parameters for describing pesticide movement through soils.  EPA-660/2-75-009.
National Environmental Research Center, Corvallis, Oregon.  May.  162 pp.

MODEL SCOPE;

     Steady-state soil water condition simulation of pesticide movement
through soil.  Isothermal conditions are assumed.  Simulation includes
pesticide concentrations with depth for soil solute and soil systems.
Transport analysis includes convection, diffusion, mechanical dispersion and
adsorption-desorption.  Output primarily consists of solute and sorbed
pesticide concentrations with time and soil depth.  A simple input dump
is included.

INPUT REQUIREMENTS;

     Control parameters provide the user with flexibility in describing the
initial state of the system.  For each run to be made, the user must specify
whether solute concentrations versus soil depth or time are to be output.  If
output is based on time, the number of time increments must also be specified.
Then a parameter is used to either generate a chemical concentration at the
soil surface or describe an initially distributed concentration.  And finally,
an initial distribution of relative solute concentration can be read in or
the program can put a square pulse distribution in initially.  Relative
solute concentration is the ratio of actual solute concentration to initial
input concentration.

     After input of the control parameters, values for flow velocity,
dispersion coefficient, adsorption coefficient, water content, bulk density,
initial concentration of pesticide, maximum soil depth considered, time of
soil application of pesticide, exponents in adsorption-desorption equations,
length of initial pesticide distribution, time increment in calculations, and
depth increment are entered.

     Other optional input data and specific details regarding format and
order are given in the previously noted reference.

SPATIAL AND TIME SCALES;

     This simulation is a one-dimensional analysis of pesticide transport
into and through a soil profile on the scale of crop root zone and the
unsaturated regions below.  Time scales may be set by user, but in general,
the simulation is of the water application-drainage interval.  In irrigated
systems, the total time scale would approximate the irrigation interval with
user supplied increments within this time frame.
                                      79

-------
COMPUTER CODE STRUCTURE;

     Program consists of a main program which coordinates data input and
output, sets time and space increments, and makes summary calculations.  Two
subroutines are included which calculate solute concentrations.  Interprogram
data transfer is accomplished with call and common statements.  Program is
written in Fortran IV and is well documented with comment statements.  Core
and execution time requirements are not given, but appear to be in the range
of 50-70 k bytes and 5-15 CP seconds per analysis, respectively.  Program
does not require tape or disk input-output devices and should be easily
adapted to computer systems other than the IBM 360 initially utilized.

BASIC MATHEMATICAL APPROACH;

     This analysis assumes a steady-state volumetric soil-water flux as
determined by the Darcy Equation.  Flow velocity is determined by dividing
the flux by the volumetric soil-water fraction.  Pesticides enter the problem
either as concentrations in applied water or as initial concentrations
distributed through the soil profile and transported through miscible dis-
placement.  The simulation of adsorption-desorption assumes the process can
be represented by a single-valued function for adsorption and desorption and
that a linear relationship exists between the adsorbed and solute pesticide
concentrations  (Freundlich equation).

     The transport of pesticide in the soil system occurs through convection,
diffusion, mechanical dispersion, and adsorption-desorption.  The transient
miscible displacement is simplified by assuming that the dispersion coeffi-
cient is independent of depth and that the Freundlich adsorption model is
adequate.  The resulting differential equation is solved by an explicit
finite difference procedure which is corrected for numerical dispersion.
Specific details of the correction procedure are found in the reference.

GENERAL;

     This model along with descriptive material is given in the previously
referenced report.  Both laboratory and field data were collected to
calibrate and verify this model.
                                     80

-------
MODEL:                                D-4

REFERENCES;

Davidson, J.M., D.A. Graetz, P.S.C. Rao, and H.M. Selim, 1977.  Simulation of
nitrogen movement, transformations, and uptake in the plant root zone.   Draft
Report to the Environmental Protection Agency, Office of Research and
Development, Athens, Georgia.  139 pp.

MODEL SCOPE;

     Simulation is of the nitrogen cycle in a crop root zone using finite
difference approximations to the flow of moisture and convective dispersion
and transport of NO^ and NH4.  Plant uptake and microbiological transforma-
tions as well as NH4 adsorption-desorption are also included.  Microbial
transformations simulated are first-order kinetic approximations of
nitrification, denitrification, mineralization and immobilization.

INPUT REQUIREMENTS;

     The user must first supply a number of control parameters to set the
number of points in the moisture content versus head relationship, initial
depth and time increments, plotter control, and the number of points various
initial conditions are specified.  Then the soil-water hydraulic conditions
are input along with initial soil nitrogen conditions, nitrogen transforma-
tion rate coefficients (nitrification, immobilization, mineralization,
denitrification}.  The Freundlich distribution coefficient along with various
soil properties such as bulk density and depth are also read in as input.
The final set of data is the description of the irrigation or rainfall event.
These parameters include duration of event, chemistry of applied water, and
evaporation.

SPATIAL AND TIME SCALES;

     This program simulates the fate of various nitrogen species in the
unsaturated crop root zone in a one-dimensional fashion.  Each model run
simulates an irrigation or rainfall event although a series can be set up
to model a longer period.

COMPUTER CODE STRUCTURE;

     The computer program consists of a main program with an input section
and seventeen subroutines and functions supplied data with call and common
statements.  Core and execution time requirements are not given but are not
expected  to restrain use on most computer systems.  The program is written
in Fortran IV and should be widely adaptable.

BASIC MATHEMATICAL APPROACH:

     The moisture flow system is simulated by solving the nonlinear partial
differential equation governing one-dimensional unsaturated flow by an
explicit-implicit finite difference approximation.  The nitrogen

                                     81

-------
transformations involve:   (1) nitrification of NH* to NOv (2) mineralization
                  T                              " J.
of organic N to NH^  (3) immobilization of both NHj and NO^ to organic N, and
(4) denitrification of NO-j to the gaseous form.  Each of these transforma-
tions was  assumed to be of the first-order kinetic type which assumes
optimal pH and temperature conditions.  Ammonium ion exchange was considered
to be an instantaneous reaction described by a linear Freundlich adsorption
relationship.  The nitrogen transformation and transport equations are also
nonlinear partial differential equations solved by finite difference
approximations.

     A nitrogen plant uptake phase was also included in the model.  The
Penman analysis of potential evapotranspiration was used in connection with a
Molz-Remson extraction pattern model to determine the ratio and distribution
of uptake.  A sink term approach was used in the uptake model and it was
assumed that root capacity for adsorption was uniform over the length of the
system.

GENERAL;

     The computer code is published in the reference cited above but is not
documented within the code itself.  Calibration and verification results are
reported.
                                     82

-------
MODEL:                                D-5
REFERENCES:
DeWit, C. T. and H. Van Keulen.  1975.  Simulation of transport processes in
soils.  Centre for Agricultural Publishing and Documentation.  Wageningen,
The Netherlands.  100 pp.

MODEL SCOPE:
     The movement of heat, salts and ions, and water in unsaturated soil is
described by a series of CSMP coded computer routines.  Transport of these
materials and energy encompasses diffusion, dispersion, and mass flow.  The
soil profile is divided into a series of compartments in which a sequential
materials balance is computed.  A number of the processes are not mutually
exclusive, and the writers do not attempt to present a clear description
for integrating the system together.  Output includes the time and spatial
distribution of the parameters.  Plotting functions are also described.

INPUT REQUIREMENTS:

     1.  Soil heat flow:  surface temperature fluctuations with time,
thickness and number of soil depth increments, thermal conductivity of the
soil compartments, and volumetric heat capacities; 2.  Salt transport:
initial moisture contents, labyrinth factors, initial solute concentrations,
inflow volume and concentration, compartment thicknesses, a dispersion
factor, and water diffusion coefficients; 3.  Ionic diffusion:  compartment
thickness, diffusion distance, total depth, water contents, labyrinth factor,
valency and initial concentrations of the ion under study, and diffusion
coefficients; 4.  Ionic transport:  initial moisture contents, labyrinth
factor, compartment geometry, dispersion factor, exchange capacities, water
inflow, diffusion coefficients, solution-adsorption site equilibrium constant,
and initial distribution of ions in the soil profile; 5.  Infiltration:
initial moisture contents, compartment geometries, initial water content at
the soil surface, hydraulic conductivity and diffusivity functions of water
contents.

SPATIAL AND TIME SCALES:
     Calculation proceeds in time increments on the order of five to ten
minutes and continues through periods of days.  The analysis is limited to
the 1-3 meter unsaturated zone and is one-dimensional.

COMPUTER CODE STRUCTURE:
     Code is written in CSMP with Fortran output.  Program is dispersed in
the reference and not presented in a single listing.  In addition, the five
segments noted are not described as an integrated model, although, as will be
pointed out next, the basic model structure is the same for each segment.
Neither core nor execution times are available.
                                     83

-------
BASIC MATHEMATICAL APPROACH;

     The soil depth is divided into a number of layers or compartments.  With
this geometry established and the initial conditions defined, the approach is
to compute the movements from one compartment to the next based on the
gradient between layers, and diffusion and dispersion functions.

     In the heat flow submodel, the heat flow rate equals the temperature
gradient between compartments multiplied by the thermal conductivity of the
soil.  The heat energy in a layer equals the temperature multiplied by the
volumetric heat capacity multiplied by the compartment size.  In order to
calculate the heat flow with time, a time distributed surface temperature is
defined along with the initial thermal characteristics of each compartment.
The dynamic nature of the system does not allow for a direct equilibrium
solution describing the heat status at each level in the soil.  The procedure
utilized is a centralized, forward integration procedure (explicit).  At any
instant in time the heat contents of all compartments are initially given.
From this information the flow rates into and out of each compartment are
computed and the volumetric heat contents are iteratively reviewed.  The
analysis can be for layered soils if their thermal coefficients can be
defined.

     The salt transport submodel operates in the same manner as the heat flow
submodel except that flow is input at the surface (rate) with salt movement
between compartments based on concentration differences (diffusion) and mass
transport with the flow.  Mass balance between compartments is maintained by
recomputing salt accumulations.  Dispersion due to tortuosity is also related
with a coefficient to the concentration gradients.

     Ionic diffusion in the absence of moisture flow is also based on the
preceding approach with gradients based on ionic concentrations.  The ionic
transport submodel includes not only ionic diffusion, dispersion, and mass
flow but cation exchange as well.  Again, the ionic composition of each
compartment is iteratively determined in conjunction with the net inflow or
outflow as was described for the heat flow submodel.  The added complexity of
the monovalent-divalent ionic adsorption processes is handled by relating the
adsorbed ions proportionally to the moisture content and exchange capacity
and then distributing the relative composition of the adsorbed ion in an
equilibrium constant approach.

     In the heat and solute movement submodels, the conductivity and
diffusion-dispersion coefficients were not dependent on the relative
concentration of the parameter being simulated.  In the case of infiltration,
however, diffusivity and hydraulic conductivity are moisture dependent.  To
simulate the moisture movement between compartments, an "average" moisture
content is determined for the gradient information.  In this submodel the
average is actually a weighted average to take into account the nonlinear
diffusivity-hydraulic conductivity versus moisture content relationships.
The flow between layers is thus the average diffusivity times the moisture
content gradient plus the average hydraulic conductivity.   The overall
procedure is to first establish the initial and boundary (inflow) conditions.
Then diffusivity and conductivity are calculated for each layer yielding the

                                     84

-------
interlayer moisture movement.  Net moisture content change and the input
water are routed through the compartments by successively computing-adjusting
until a stable mass balance is achieved or until changes from iteration to
iteration are small.

GENERAL:
     Although the model is detailed in the reference, a potential user is
still left with a considerable effort to integrate the parts into a useable
program.  Fortran users will be able to develop such models with about as
much effort.
                                      85

-------
MODEL;                                D-6

REFERENCES;

Frissel, M. J. and P. Reiniger.  1974.  Simulation of accumulation and
leaching in soils.  Centre for Agricultural Publishing and Documentation,
Wageningen, The Netherlands,  pp. 70-84.

MODEL SCOPE:

     Simulation is of diffusion, dispersion, adsorption-desorption, mass
transport, and volatilization of herbicides in a homogeneous soil profile.
Output includes the time and spatial distribution of the herbicides in the
soil system.

INPUT REQUIREMENTS;

     Two options are available:  (1) instantaneous equilibrium  (no
adsorption); and (2) noninstantaneous equilibrium (adsorption-desorption).
Both model options require volumetric moisture contents, bulk density, gas
and liquid phase fractions, distribution ratios (soil/water, water/gas), soil
profile depth, diffusion coefficient for gaseous phase, surface flux and
concentration, and first-order decomposition rate coefficients for the
herbicides.  The noninstantaneous option also requires adsorption-desorption
rates and water system diffusion rates.

SPATIAL AND TIME SCALES:
     Simulation is of one-dimensional analysis of a one meter unsaturated
soil profile.  Calculations are presented at intervals specified by the user
in hours or days.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The flux of the herbicides due to water movement is described basically
by the compartment modeling approach noted in other CSMP models.  For the
instantaneous equilibrium model, herbicides migrate in both the water and
gaseous phases (only diffusion flux for the gaseous phase).  Decomposition of
the herbicide is a first-order reaction based on a decay rate.  The analysis
describes the convection and diffusion flux in the water phase and the
diffusion flux in the gaseous phase.  In the noninstantaneous equilibrium
model, a linear adsorption reaction is assumed and added to the instantaneous
equilibrium approach.
                                     86

-------
GENERAL;

     The mathematics for the herbicide models are well defined in the
reference.  The models are listed in the reference.  Field verification is
not reported.
                                     87

-------
MODEL;                                D-7

REFERENCES;

Frissel, M. J. and P. Reiniger.  1974.  Simulation of accumulation and
leaching in soils.  Centre for Agricultural Publishing and Documentation,
Wageningen, The Netherlands,  pp. 54-69.

MODEL SCOPE:
     Simulation of diffusion, dispersion, and mass transport of Na , K ,
Ca  , and Mg++ in soils.  Model considers nonlinear adsorption-desorption.
Output includes the time and spatial distribution of the salinity cations in
the soil solute, adsorbed phase, and in deep percolation.  The soil is
considered homogeneous.

INPUT REQUIREMENTS;

     Input consists of depth distributed and exchangeable solution Na , K ,
Mg  , Ca++, moisture content, Na+/K+, Ca++/Mg"H", and adsorbed cation fraction/
cation exchange capacity ratios, inflow rates and concentrations, fixed K
and fixation and release rates.

SPATIAL AND TIME SCALES;

     Spatial scale is a one-dimensional analysis of a one meter or so depth
of unsaturated soil.  Calculations are reported in hours.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The diffusion, dispersion, and mass transport are simulated by the CSMP
compartmental model described previously.  In addition, the cation exchange
and potassium fixation are included to adjust the transport results.  The
exchange reaction is assumed to occur instantaneously.  A Gapon equation,
along with nonlinear functions (Vanselow), are incorporated, but the principal
approach is one developed by the authors.  Potassium fixation is a first-
order reaction depending on relative concentration in solution and on the
soil.

GENERAL:
     More detail concerning the exchange processes is left to the interested
reader.  The input can be time distributed, but it is not clear how the
moisture is redistributed since moisture content-conductivity relationships
are not included.  Program is listed in the cited reference.
                                      88

-------
MODEL:

REFERENCES;

Frissel, M. J. and P. Reiniger.  1974.  Simulation of accumulation and
leaching in soils.  Centre for Agricultural Publishing and Documentation,
Wageningen, The Netherlands,  pp. 12-20.

MODEL SCOPE:
     Simulation of diffusion, dispersion, and mass transport of completely
soluable compounds in the soil solute.  Output includes the time and spatial
distribution of the solute concentration and deep percolation.  The soil is
considered homogeneous.

INPUT REQUIREMENTS;

     Input to the model includes depth of soil profile, number of compartmental
soil layers, surface flux and inflow concentration, diffusion constant, anion
exclusion ratio, initial moisture content, tortuosity and dispersion coeffi-
cients with soil depth, output frequency, and initial distribution of solutes.

SPATIAL AND TIME SCALE:
     The spatial scale is an approximately one meter, one-dimensional view of
an unsaturated  soil profile.  Calculation times are on the order of minutes
with output printed at selected aggregate intervals of hours.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The soil profile is divided into 40 compartments or layers with initial
conditions defined in each layer.  The diffusion coefficient for the boundary
between two  layers is computed using mean values of moisture content, tor-
tuosity, and dispersion values.  Diffusion rates are based on Pick's Law
which calculates an average concentration gradient between layers and multi-
plies it by  the diffusion coefficient to arrive at the flow.  Mass flow
is also based on average concentration and moisture flow rates.  In this
model, water flux is based on a steady-state infiltration.

     The basic mathematical approach involves first calculating the movement
rates between compartments.  The rate of change in compartmental solute
concentration is then computed and integrated using a standard CSMP semi-
parallel method such as a fourth order Runge-Kutta method.  The rates are
recomputed and the iteration repeated until the compartment concentrations
have stabilized.
                                     89

-------
GENERAL;

     The program is set up for tritiated water.  For other anions (Cl ,
NO3), the rates are multiplied by the exclusion ratio.  Code is presented in
the cited reference.  Verification is not reported.
                                     90

-------
MODEL;                                D-9

REFERENCES;

Hagin, J. and A. Amberger.  1974.  Contribution of fertilizers and manures to
the N- and P-load of waters, a computer simulation model.  Report to the
Deutsche Forschung Gemeinschaff.  January.  123 pp.

MODEL SCOPE;

     This detailed simulation model addresses the conditions which induce
nitrate and phosphate pollution of receiving waters by fertilizers and animal
wastes.  The model evaluates the short term chemical, microbiological, and
physical processes in the soil and then aggregates these results in the
prediction of longer term outflow of nitrates and phosphates from the system.
Nitrogen is assumed to originate as NO^, NH^, and organic material having
differing C:N ratios.  Simulation of subsequent nitrogen transformations
include mineralization of organic-N to NH^, nitrification of NH| to NO^
(assumed to occur in top 10 cm of the soil), and denitrification.  Also
included are the transport and plant uptake of NO .  Phosphate is added to
the soil surface and then moves through soil erosion.  Model does not have
either input dump or error flag systems.  Input is accomplished by internally
punched statements.

     The nitrogen model is divided into the following sections:  (1) water
movement; (2) heat movement; (3) ammonification of organic-N;  (4) oxidation
of ammonium; (5) oxygen movement;  (6) denitrification; and (7) nitrate
movement.

     The nitrogen model output includes hard copy plotting routines to
present water content versus depth at different times, nitrate profiles,
surface temperature profiles, cumulative plant uptake of nitrate, oxygen
profiles, and total nitrate present.

     The phosphate model consists of a section on water flow and one for soil
detachment,  transport, and deposition.  Graphical output includes phosphorous
erosion and sedimentation along field slope or outflow.

INPUT REQUIREMENTS;

     The nitrogen model consists of 7 sections.  Section 1, Water Movement,
input includes rainfall and irrigation duration and intensities, class A pan
evaporation, crop coefficients, soil hydraulic conductivity versus moisture
content, conductivity at soil-atmosphere interface and various soil segments,
diffusivity in various soil segments, depths and number of homogeneous soil
layers, initial soil moisture conditions,  and various root distribution and
strength data.

     Section 2,  Heat Movement Data, includes air temperature, heat transfer
coefficients, capacities,  and initial conditions, and temperature of
irrigation or rain water.


                                     91

-------
     Section 3, Ammonification of organic nitrogen, requires the rate versus
pH, temperature, and moisture dependent coefficients, percent organic matter
added to soil as cellulose and other forms, immobilization-mineralization
rates for Lignin, hemicellulose, sucrose, and cellulose, and additions or
initial amounts of NO -N.

     Section 4, Oxidation of Ammonium, input involves applied NH^, relative
mineral composition of soil and soil surface layers, soil material density,
initial NO- distribution, available nitrification bacteria, surface layer
nitrification rate constants and coefficients, soil pH and temperature.

     Section 5, Oxygen Movement, data include aerobic respiration rates,
initial oxygen content in water throughout profile, and various other
coefficients such as diffusion relationships.

     Section 6, Denitrification, input includes carbon percentage of existing
and added materials, denitrification rate coefficients and constants, and
various initial conditions.

     Section 7, Nitrate Movement, analysis requires input from previous
sections, plus initial nitrate distributions, and diffusion and dispersion
factors.  Additional tortuosity and root uptake coefficients are also
required.

     These are by no means an exhaustive list and some misinterpretation of
the code may result in listing an input where the actual data is computed.

     For the phosphate part of the model, the water movement section
requires geometries of soil subdivisions, rain and irrigation durations and
intensities, surface hydraulic parameters, soil infiltration rates, porosities,
and erosive parameters, slope of surface, etc.  The second section requires
data describing added nitrogen and phosphorus.

TIME AND SPATIAL SCALES;

     Simulation is of a single or representation site on an irrigated or rain
fed field.  Time increments vary somewhat in the model from an almost con-
tinuous system for nitrogen transformation to daily events for irrigation or
precipitation to weekly periods for parts of the evapotranspiration reactions.
In general, the time scale is approximately a daily breakdown of the irriga-
tion season.  Model is one-dimensional.

COMPUTER CODE STRUCTURE;

     CSMP

BASIC MATHEMATICAL APPROACH;

     (1)  Water Movement.  Water is added to the upper soil layer through
either daily inputs of precipitation or irrigation.  Transpiration is
determined with class A pan data multiplied by plant growth coefficients.
The uptake by the roots represented by the transpiration is based on the

                                      92

-------
assumption that the uptake is proportional to the transpiration rates and
distributed between soil layers according to relative root activity.  A limit
of 20 atmospheres is set on the moisture withdrawal in a layer.  Water loss
through soil surface evaporation is simulated for flooding, no-flooding,
irrigation, and no-irrigation.  Under no-flooding, no-irrigation, the evapora-
tion is actually the difference between evapotranspiration as computed using
pan data and transpiration.  If the upward flow of water by diffusion is less
than the previous difference, the evaporation is set at the smaller value.
In the flooding case, the evaporation is equated to the class A pan values
but for irrigation the corresponding value is assumed to be some minimal
value.

      (2)  Heat Movement.  The temperature of each soil layer is predicted by
dividing the volumetric heat content by the volumetric heat capacity.  Heat
into and out of each layer due to moisture movement and transpiration is
determined by multiplying flow rates by temperature.  Heat conduction between
layers is determined by multiplying the average heat conductivity of the two
adjacent layers by the temperature gradient between layers.

      (3)  Ammonification of Organic-N.  Mineralization and immobilization of
nitrogen is assumed to occur in the upper 10 cm of the soil in a first-
order reaction in which the mineralization rate is proportional to the
available organic N.  A distinction is made in organic materials according  to
their most common carbohydrate, lignin, cellulose, hemicellulose, and
sucrose.  The C:N ratios delineating the two processes vary with the
carbohydrate.  It should also be noted that additional reactions such as
humification, the pH and moisture effects, and lignin inhibition are also
considered.

      (4)  Oxidation of Ammonium.  The nitrification rate equation used in
this model is a first-order sigmoid curve.  Only soil pH is evaluated
directly in the model; temperature, moisture, and oxygen concentrations are
assumed at their optimum values.  The influence for each of these parameters
is derived from literature reports.

      (5)  Oxygen Movement.  Oxygen diffusion is modeled with Pick's second
law by assuming simultaneous diffusion in both air and water in the soil.  It
is also assumed that instantaneous equilibrium between gaseous and dissolved
oxygen occurs at any soil depth.  Oxygen utilized in the oxidation of ammonium
and organic carbon is also evaluated.  This part of the model is one of the
more complex included in the simulation.  The reader is directed to the
reference for specific details regarding definition of the diffusion coeffi-
cients, interlayer transfer, and consumption by organic carbon and ammonium
processes.

      (6)  Denitrification.  Probably the least developed of the nitrogen
transformation simulations is denitrification.  Hagin and Amberger  (1974)
treat the process as a first-order reaction with a rate coefficient
determined by multiplying coefficients for temperature, pH, available organic
matter, and oxygen-moisture.
                                     93

-------
     (7)  Nitrate Transport and Uptake.  Nitrate is assumed to occur in the
profile as an inflow with irrigation water, nitrification in the first soil
layer, and as applied fertilizers.  Nitrate losses consist of deep percola-
tion, plant uptake and denitrification.  Interlayer migration occurs either
by diffusion or mass flow.  Diffusion is proportional to the concentration
gradient between layers and inversely proportional to the distance between
layers.  Mass flow occurs with the moisture movement but is adjusted with a
dispersion coefficient.  The uptake of nitrate is calculated by multiplying
the volume of moisture uptake by the soil water nitrate concentration  (with
unit conversions) but limited to some upper limit.  Finally, nitrates can
accumulate in a soil layer or be leached depending on a mass balance
computation.

     For the transport of phosphates in surface runoff and erosion, the upper
soil layer is assumed sufficiently permeable that runoff would not occur
until it is saturated.  The depth of surface runoff is computed much like the
surface irrigation advance analyses, i.e., a mass balance between the surface
and infiltrated moisture.  The erosion mechanism includes soil detachment by
rainfall and runoff and transport capacities of both rainfall and runoff.
The individual processes utilized in this phase of the model are presented in
the technical literature.

GENERAL;

     This is one of the more detailed treatments of the nitrogen-phosphate
system in irrigated agriculture, but most of the model is developed from data
published in the literature.  The calculations presented are likewise based
on reported data.  Consequently, the "general" nature of the model may remain
somewhat in question until it is applied and verified against field data
taken at one location.  It should be useful in evaluating the interrelation-
ship in this complex system and is therefore primarily a research tool.  The
code is included in the cited reference.
                                      94

-------
MODEL;                                D-10

REFERENCES;

Hillel, D.  1977.  Computer Simulation of Soil Water Dynamics.  International
Development Research Centre, Ottawa, Canada,  pp. 131-134.

MODEL SCOPE:
     Simulation is of the water budget and noninteracting solute transport in
the unsaturated crop root zone.  Model includes root extraction, root growth,
and the spatial-time distributed flux of water and solute.

INPUT REQUIREMENTS;

     Parameters are specified defining the geometry of the soil and root
system, hydraulic properties of soil and roots, daily transpiration demand,
hydrodynamic dispersion coefficient, and water potential at the plant crown.
Tables of relative root length in the vertical profile, initial moisture and
solute conditions, suction, wetness and conductivity relationships, and
wetness-tortuosity relationships are specified.

TIME AND SPATIAL SCALES:
     This model is a one-dimensional treatment of the unsaturated root zone
determined for one second intervals aggregated to hourly and daily summaries.

COMPUTER CODE STRUCTURE:
     CSMP

BASIC MATHEMATICAL APPROACH;

     The compartmentalized soil system is assumed to be in a state similar to
post-irrigation or post-rainfall by specifying the initial distribution of
soil moisture.  Moisture movement between layers is based on calculating an
average hydraulic gradient and conductivity between soil layers or compart-
ments.  Then intercompartmental flux is used iteratively to recompute
moisture content.

     The uptake segment of the model involving developed or growing root
systems is based on distributing the transpiration through the soil profile.
The relative extraction in a compartment is the ratio of the difference
between the moisture potential in the compartment and the crown potential to
the sum of the hydraulic resistance of the soil  (conductivity x length of
flow path) and the compartment's root resistance.  The crown potential is
identified by successive approximation.  If root birth and extension are
considered, exponential functions are used to simulate them so the rooting
pattern can be updated at each iteration.
                                      95

-------
GENERAL;

     Verification trials are not reported although results are presented
illustrating the model results.  The code is listed in the reference.
                                      96

-------
MODEL;                                D-ll

REFERENCES;

Margheim, G.A., 1967.  Predicting the quality of irrigation return flows.
Unpublished M.S. Thesis, Civil Engineering Department, Colorado State
University, Fort Collins, Colorado, December.  57 pp.

MODEL SCOPE:

     Simulation is of the salinity constituents in flows entering the
groundwater basjn under irrigated croplands.  Parameters include the concen-
trations of Ca  , Mg  , Na , and 504 in the flows entering the water table
from deep percolation.  The computer code does not contain an input dump or
a printing of intermediate results, but does contain partial error flags for
iteration limitations.

INPUT PARAMETERS;

     No control parameters are utilized.  Input data_include the irrigation
water concentrations of Ca++, Mg++, Na+, Cl~, and 304, and the concentrations
of exchangeable Na+, Ca++, and Mg++ in the soil.  Input also requires concen-
trations of soluable Mg"1"1", Na+, Ca++, and 304 in the soil, solubility product
of gypsum, and cation exchange constants for homovalent and monovalent-
divalent exchange equations.

SPATIAL AND TIME SCALES;

     Analyses are at a single or representative site in an irrigated field
(accomplished by assuming uniform irrigation water applications).  Steady-
state moisture flux is assumed so time scale is set by user.  However, an
irrigation interval is generally selected.

COMPUTER CODE STRUCTURE;

     Internal program linkage is facilitated by subroutines with data and
results transferred through call statements.  Program can be easily adapted
to other computer systems since tape and disc functions are not performed.
Core requirements are less than 40 k bytes while execution times are not
known (probably 10-15 seconds per run).  Model is programmed for CDC 6400
using Fortran IV language.

BASIC MATHEMATICAL APPROACH;

     Under assumptions that gypsum is the only slightly soluble salt present
in the soil system, the model first computes the ion concentrations in water
percolating below the root zone and then the concentration of water inter-
cepted by a drainage system.  Effects of evapotranspiration are not directly
considered in the soil chemistry phase.  Computations in the soil chemistry
phase involve the following iterative procedures controlled by iteration to
iteration changes in calcium concentrations.  A volume of irrigation water
equal to the stored soil moisture volume is equilibrated with the soil-CaS04

                                     97

-------
system using the Gapon equation for exchange processes, and the
Debye-Huckel gypsum solubility relations for solution equilibrium.  The
iterative procedure is repeated once the successive approximation procedure
has converged for other volumes of irrigation water until the desired leaching
is accomplished.  Then, using a drainage equation developed by R.E. Glover,
the volume of return flow is computed.  The model assumes groundwater quality
is poorer than deep percolation and mixing does not occur.  An interface is
thereby formulated from which the relative displacement is determined to
arrive at the outflow mixture.

GENERAL:

     Code is available in reference cited above.  However, code is not
documented nor verified against field data.
                                     98

-------
MODEL:                                 D-12

REFERENCES:

Childs,  S. W.f  and R.  S.  Hanks.   1975.  Model of  soil  salinity effects on
crop growth.   Soil Sci. Soc. Amer. Proc., 39:617-622.

Nimah, M. D.,  and R. J. Hanks.   1973.  Model for  estimating soil water,
plant, and atmospheric interrelations:  I. Description and Sensitivity, and
II. Field test of model.   Soil Sci. Soc. Amer. Proc.,  37:522-527 and
37:528-532.

Melamed, J.  D.,  R.  J.  Hanks, and L. S. Willardson.  1977.  Model of salt flow
in soil  with a source-sink term.  Soil Sci. Soc.  Amer. Proc., 41:29-33.

MODEL SCOPE;

     The model described  here is actually the result of several individual
models.  Nimah and Hanks  (1973)  added a plant root extraction system to an
earlier  one-dimensional moisture flow model.  Childs and Hanks (1975)
modified the Nimah and Hanks  (1973) model for crop growth, salt flow, and
irrigation uniformity.  Then, the "Dutt Model" described herein by Dutt et
al. (1972) and a simple source-sink model to account for chemical reactions
were added.   The simpler  salt system was added to the  water model by Melamed
et al.  (1977).   Because the more complex Dutt model is described elsewhere,
the version  will not be detailed here.

     The model  simulates  the one-dimensional flow of water in the transient,
unsaturated  soil profile  including the spatial distribution of root uptake in
response to  surface evapotranspiration demand.  Salts  added to the system in
the irrigation  water or that exist in the solid form in the soil structure
are transported through the system to the groundwater.  below.  Dissolution and
precipitation of salts in the root zone are simulated  as a source-sink
process  rather  than by chemical  equilibrium.  Ionic exchange and biological
transformations  are neglected.   Output describes  the spatial and time
distributed water-salt system including evapotranspiration.

INPUT REQUIREMENTS;

     The primary data  necessary  for the soil moisture  flux are the hydraulic
conductivity-water  content and pressure head-water content relations, air dry
and saturated moisture contents,  root water potential  below which wilting
occurs, depth distributed root distribution, initial moisture distribution,
potential evapotranspiration at  the surface and water  inputs (irrigation or
rainfall), osmotic potential of  the irrigation water and of soil (depth
distributed), the presence of a  lower boundary water table, and crop cover.
The salinity simulation requires  the electrical conductivity of the irriga-
tion water and solid soil  salts,  diffusion and dispersion coefficients, and
initial solute distribution (electrical conductivity).
                                     99

-------
TIME AND SPATIAL SCALES:

     The model is a one-dimensional analysis of the root zone.  Time scales
vary in length from an irrigation interval to several years.

COMPUTER CODE STRUCTURE;

     Programs are written in Fortran IV with common and call data transfer to
subroutines.  Execution time depends on the length of the real time simulation.
Core requirements appear moderate (50-70 k bytes).

BASIC MATHEMATICAL APPROACH;

     The model involves a moisture flow simulation, a salt transport system,
and a simulation of the dissolution-precipitation processes  (source-sink
term).  The moisture flow equation is the moisture content form of the
unsteady flow equation developed by writing Darcy's Law and mass conservation
equations in one dimension.  This primary flow equation is modified by the
inclusion of a root extraction term and solved with a implicit finite
difference technique.  The essential features of the root extraction term are
described by Nimah and Hanks (1973).  (See also the model of Feddes and
Zaraday, 1977.)   Flow to the roots is based on the difference between the
root potential plus root resistance and the combined soil water and osmotic
potentials.  This value is then multiplied by the hydraulic conductivity and
a depth distributed rooting pattern.  The root potential is determined by
maximizing transpiration as was the case in the model by Feddes and Zaraday
(1977)  (which was based on this model).

     The solute transport relationships are based on the equation of Bresler
(.1973) which includes diffusion and hydrodynamic dispersion.  This equation
is solved using a tridiagonal matrix solution of a Crank-Nicholson finite
difference method.

     The source-sink term accounting for precipitation and dissolution of
salts is an on-off linear equation based on the difference the salinity
concentration of the soil solution and a value at which precipitation will
occur.  The equilibrium concentration must be evaluated in the field.

GENERAL;

     The model which has been calibrated against field data shows good
predictive agreement with observed field data once calibrated.  The code
should be requested from its authors.
                                     100

-------
MODEL;                                D-13

REFERENCES:

Smajstrla, A.G., D.L. Reddell, and E.A. Hiler.  1974.  Simulation of miscible
displacement in soils.  Paper 74-2015 presented at the Annual Meeting of
ASAE, Stillwater, Oklahoma.  June.  31 pp.

MODEL SCOPE:
     Simulation is of vertical miscible displacement of noninteracting
solutes in unsaturated, homogeneous, isotropic soils.  Analysis includes the
effects of convection, ionic diffusion, and mechanical dispersion.  Output
describes the time and depth distributed values of water content and solute
concentrations during an infiltration event.

INPUT REQUIREMENTS;

     The first group of input data includes the number of soil layers being
considered, number of moving points, a control for infiltration versus
redistribution phases, number of points in the capillary pressure versus
water content table and the number of points in the water content-hydraulic
conductivity table.  Then the pressure versus water content and water content
versus hydraulic conductivity tables are read in.  Next, the grid character-
istics are defined  (initial moisture and solute state, thickness, and loca-
tion of moving points).  Following these data, the saturation moisture
content, saturated hydraulic conductivity, surface pressure during infiltra-
tion, total intake of water and salt, salt concentrations and tracer depths
allowed, simulation time and time infiltration stops, output intervals, and
several increment values not specified are assigned.

SPATIAL AND TIME SCALES:
     Model evaluates the one-dimensional transport of soil solutes in the
unsaturated zone, presumably in an irrigated situation.  The overall time
scale is the infiltration-redistribution period, something on the order of
one to five days.  Time intervals on the order of tenths of hours are used to
solve the flow equations.

COMPUTER CODE STRUCTURE:
     The computer code is written in Fortran IV for IBM 360/65 computer but
appears readily adaptable to most other kinds as well.  A main program serves
only to control the order of calculations by sequentially addressing various
subroutines.  Data is maintained in common blocks.  Core storage requirements
do not appear to exceed 50-70 k bytes but execution time would be highly
variable depending on the length of simulation and resolution in the grid
system.
                                     101

-------
BASIC MATHEMATICAL APPROACH;

     This analysis assumes a homogeneous, isotropic soil with respect to both
physical and hydraulic properties.  In addition, it is assumed that head
versus moisture content and hydraulic conductivity are single valued
functions, ions are noninteracting and conservative, and the viscosity of the
water is constant.  The nonsteady Darcy equation and the convective-dispersive
solute transport equations are solved simultaneously.  The evaluation of both
convective and dispersive processes is made using the method of character-
istics which effectively separates the two processes and thereby eliminates
numerical dispersion.  Solution of the moisture flow equation is accomplished
with an explicit finite difference technique, as was the solution for dis-
persion effects.  A system of moving points within the grid was used with the
method of characteristics for calculating the moving position of the solute.
A more detailed, description of the numerical procedure can be abstracted
from the reference cited above.

GENERAL;

     The program is well written and documented.  Although the analysis would
be limited to the conservative ions like chloride, the results would be
useful in predicting selected aspects of the irrigation return flow system,
for example, leaching fractions.  Code is listed in the cited reference.
                                      102

-------
                                APPENDIX E

                DESCRIPTION OF WATERSHED HYDROLOGY MODELS
MODEL:                                E-l
REFERENCES:
Hillel, D.  1977.  Computer Simulation of Soil-Water Dynamics.  International
Development Research Centre, Ottawa, Canada,  pp. 155-198.

MODEL SCOPE:
     Simulation is of the hydrology of a sloping agricultural field.
Processes modeled include precipitation, infiltration, runoff, evaporation,
redistribution, deep percolation, capillary rise from the water table, and
groundwater drainage.  Plant uptake and transpiration are not included.
Output depicts the time and spatial distribution of water in the system.

INPUT REQUIREMENTS;

     The geometry of the system is defined by stating the number of columns,
number of layers per column, width of each column, and vertical dimensioning.
Initial conditions such as saturated hydraulic conductivity and initial
moisture distributions are detailed.  The runoff hydraulic parameters are
input, including roughness factors for the Manning's equation and surface
slope.  Boundary conditions are defined to depict the potential evaporation
rate, water table elevation at its outflow to the drain and rainfall inputs.
As in other unsaturated soil models, the relationships between water content,
hydraulic conductivity, and suction must be known.

TIME AND SPATIAL SCALES;

     The scale is a two-dimensional evaluation of a field sized system
extending vertically to an impervious layer.  The principal event is the
period immediately preceding and shortly after a rainfall (until equilibrium
conditions occur).

COMPUTER CODE STRUCTURE;

     CSMP
                                     103

-------
BASIC MATHEMATICAL APPROACH;

     Water added to the soil surface is divided into infiltration, surface
storage and runoff.  Infiltration is calculated by dividing the vertical
unsaturated soil profile into layers and then routing flows through it using
an iterative computation of contents, gradients, rates, and new contents.
Runoff is handled by dividing the field into a number of vertical columns and
then routing the surface water down slope with a successive kinematic wave
equation.  At the bottom of the slope the time distributed runoff is yielded
by the last computations.  In each column, the infiltration crossing the last
compartment in the unsaturated zone is the inflow to the groundwater basin.
A water table height increase is distributed through each column and then the
saturated soil properties along with hydraulic gradients are used to
calculate a lateral groundwater flow for the bottom of the slope.

GENERAL:
     No field verification is presented.  Code is listed in cited reference.
                                     104

-------
MODEL:                                E-2
REFERENCES:
Makkink, G. F. and H. D. J. Van Heemst.  1975.  Simulation of the water
balance of arable land and pastures.  Centre for Agricultural Publishing and
Documentation, Wageningen, The Netherlands.  79 pp.

MODEL SCOPE:
     Simulation is of water movement and storage in a cropped field and the
saturated-unsaturated soil beneath.  The processes modeled include snow
precipitation and melt, canopy interception and evaporation, soil surface
evaporation and infiltration, unsaturated and saturated flow, evapotranspira-
tion, micellar flow and plant withdrawal and hydration-dehydration.  Output
involves the time and spatial distribution of water in the system.

INPUT REQUIREMENTS;

     The data listed as required input include the date of the last balance
period and the current date, global and extraterrestrial solar radiation, a
"place factor", air and dew point temperatures, relative humidity, wind
speed, crop height, and depth of root zone, all for each day of the analysis.
However, a number of data not listed as input are required nevertheless.  For
instance, the conductivity and moisture versus suction relationship, initial
conditions, infiltration, and soil and crop characteristics are internal to
the program as DATA statements.

SPATIAL AND TIME SCALES:
     The analysis is of a large crop area, but the evaluation is made at one
vertical "average" location.  Time increments of 0.2 days are used in the
calculation with output based on daily intervals.  Model is one-dimensional
and steady-state.

COMPUTER CODE STRUCTURE:
     Code is written in Fortran IV for the CDC series computers.  A main
program with common linked subroutines and functions comprise the basic code
structure.  Although neither core nor time requirements are discussed, they
are obviously small enough for most computers.

BASIC MATHEMATICAL APPROACH:
     The plant-soil system is divided into "compartments" representing the
various forms of moisture storage in the system.  The computational system
involves an iterative evaluation of compartment storage, the rate of moisture
movement between compartments, and the calculation of contents (CONTENT at
time t = CONTENT at time t-1 + RATE times Delta Time).  The storage compart-
ments include snow  (existing as solid precipitation or heavy frost and
flowing in the system by evaporation and melting), adhering water which stems
from canopy interception (this moisture evaporates, drops into pools, or

                                    105

-------
drops onto and infiltrates unsaturated soil).  Pools filled with precipitation
at rates greater than the saturated hydraulic conductivity of the soil
(released from pools by evaporation and infiltration), the unsaturated soil
profile  (transmitting water to the saturated zone below or to the atmosphere
above through transpiration and evaporation), the saturated zone  (venting to
transpiration, groundwater flow, water table storage changes, and flow to the
unsaturated region through capillary use) and micellar water taken up by and
stored in the clay structure.  Subsets of the mass balance in the unsaturated
zone are developed for the upper soil profile contributing to evaporation and
the region of plant uptake.  There are also storage balances for the moisture
stored above and below 16 atm and field capacity, respectively.  The basic
rate processes are evapotranspiration and soil moisture redistribution.
Evaporation and transpiration are computed with modified forms of the Penman
equation.  Unsaturated flow occurs when the soil moisture level exceeds field
capacity.

GENERAL:
     Although this model has many aspects of a general nature, its existing
format is entirely site specific to the case studies it was validated against.
To apply this model to other situations, much more data are needed.  The code
is presented in the reference cited above.
                                     106

-------
                                APPENDIX F

              DESCRIPTION  OF WATERSHED HYDRO-QUALITY MODELS
MODEL;                                F-l

REFERENCES;

Dutt, S.R., M.S. Shaffer, and W.S. Moore, 1972.  Computer simulation model of
dynamic bio-physicochemical processes in soils.  Technical Bulletin 196,
Agricultural Experiment Station, University of Arizona, Tucson, Arizona.
October.  101 pp.
Shaffer, M.S., R.W. Ribbens, and C.W. Huntley.  1977.  Prediction of mineral
quality of irrigation return flow, Volume V.  Detailed return flow salinity
and nutrient simulation model.  EPA-600/2-77-179e.  Robert S. Kerr Environ-
mental Research Laboratory, Office of Research and Development, U.S. Environ-
mental Protection Agency, Ada/ Oklahoma.  August.  228 pp.

MODEL SCOPE;

     This model simulates chemical and physical processes associated with
agricultural lands drained by subsurface drainage systems.  Basic input data
involves the field application of irrigation or precipitation with associated
salts and nitrogen nutrients.  Model output is the prediction of drainage
quality with evaluations at intermediate points in the plant-soil-aquifer
system.  Program includes intermediate output and various error flags.

INPUT DATA;

     This is a very large computer model well documented by Shaffer, et al.
(1977).  As a consequence, the input requirements are elaborate.  The
following input data summary is given by Shaffer, et al., (1977):

      (1)  Drainage parameters:  drain spacing and depth, depth to impermeable
          barrier, envelope size, saturated hydraulic conductivity, porosity,
          and specific yield or storage coefficients;

      (2)  Unsaturated flow parameters:  the hydraulic conductivity versus
          moisture content function, and pressure head versus moisture
          content relations;
                                           -L>J-    -4-^»    -4-    —     mm    •*•"
      (3)  Soil chemistry data:  meq/1 of Ca   , Mg  , Ma , CO3, HC03, Cl ,
          304, NH4, NO3, and urea in soil extract, pH, cation exchange
          capacity, gypsum content, presence of lime, bulk density, organic


                                     107

-------
          nitrogen, carbon-nitrogen ratio, and the nitrifier population-salt
          response relationship;

     (4)  Crop information:  rooting depth and distributions,
          evapotranspiration rates as read in or calculated using an
          irrigation scheduling program, and plant uptake of nitrogen;

     (5)  Water application data:  irrigation schedules and amounts (or
          precipitation amounts and timing) and an effective precipitation
          relationship;

     (6)  Fertilizer data:  fertilization schedules, amounts of NO3, NH^,
          Urea, Ca"*"+, SO^, and CO^, application depth and organic nitrogen-
          plowed in with corresponding C-N ratio, and application ratio;

     (7)  Irrigation water analyses:  concentrations of NH4, NO3, Ca++, Na+,
          Mg++, HCO3, Cof, Cl , and 304; and

     (8)  Miscellaneous soils data:  soil temperature and coefficients
          relating to nitrification and denitrification.

SPATIAL AND TIME SCALE;

     The spatial scale is a two-dimensional evaluation of the vertical region
from the soil surface to the drainage system.  Time scales are mixed depending
on the calculation being made.  Time varies from fractions of a day to an
irrigation interval and even up to an irrigation season.  The model is
transient in nature.

COMPUTER CODE STRUCTURE;

     The primary code structure is in a main program-subroutine format
utilizing common, call, and tape file data transfer.  The size of this model
requires extensive use of overlay systems if the entire package is used.
Consequently, core requirements are only 140-150 k bytes at any instant.
Time for execution varies from one to two minutes up to several minutes
depending on the scope of the job.  Although the basic program language is
Fortran IV, this model is not easily adapted in its entirety.  Except for
some planning use by the Bureau of Reclamation, this is primarily a research
tool and not applicable by most planning groups.

BASIC MATHEMATICAL APPROACH;

     This model is comprised of individual models which extend across several
disciplines.  Consequently, there is more likely to be more interest in the
various submodels rather than the total package.  The program is linked
together with detailed overlay systems and numerous by-pass options provide
substantial flexibility.  The interested user can find excellent explanatory
comments in the references noted above.  However, to identify the model's
basic characteristics, this section will be divided by major submodels.
                                     108

-------
     Irrigation scheduling.  Unlike the irrigation scheduling program
available to schedule and update irrigations, this model predicts irrigation
schedules and amounts using historical data exclusively.  Based on geograph-
ical data, solar radiation, and temperature, the program first computes
potential evapotranspiration for alfalfa using the Jensen-Haise equation.
Adjustments are made for seven other crops using standard growth stage coeffi-
cients.  Then, a simple mass balance of the root zone is developed.  Soil
moisture is allowed to be depleted to an "allowable moisture depletion"
before an  irrigation is scheduled.  Coefficients for various irrigation
efficiency and losses produce deep percolation and surface runoff estimates.
The output from the model is recorded on cards punched for subsequent use in
the unsaturated moisture-chemistry models  (infiltration) and the drainage
calculation  (deep percolation).

     Soil moisture movement.  The unsaturated flow submodel is a finite
difference simulation of the transient moisture flow equation in one dimen-
sion.  Hydraulic conductivity is taken as a function of moisture content as
is pressure head.   Simulation includes infiltration, redistribution,
drainage, and plant uptake.  Plant uptake is determined from information of
total uptake which is distributed throughout the soil profile in proportion
to an average root distribution.

     Unsaturated soil chemistry.  The complex soil-water system is simulated
from a chemical-biological standpoint by a group of subroutines designated as
the unsaturated chemistry program.  The program has two primary components:
(1) simulation of the soil nitrogen system; and  (2) the soil salinity series.
This model is the basic programming package reported by Dutt, et al.  (1972).

     The nitrogen phase of the program considers urea hydrolysis, nitrification
of NH^-N, net mineralization-immobilization of both organic-N and NH4-N,
immobilization of N03-N, plant uptake of NH^-N and N03-N, and denitrifi-
cation.  In  addition, the NH^-N reaction in the cation exchange system is
considered.  With the exception of the cation exchange process, the nitrogen
simulation assumes a first-order kinetic approach using regression functions
for the transformation rates.  Denitrification is evaluated using zero-order
kinetics.  Nitrogen uptake may be evaluated by either of two methods.  First,
the user may specify total N uptake and root distribution which yield the
plant N uptake in each soil segment per unit of time.  In this case, the user
specifies the fractions of NO3, and NH^ uptake.  The second assumes uptake is
proportional to water uptake.  Output from these segments of the model
predicts mass distribution of the various nitrogen components.

     The salinity phase of the models evaluates ion exchange, solution-
precipitation of gypsum and lime, formation of undissociated ion pairs, and
transport of the soluble species  (Ca++, Mg"4"4", Cl  , Na+, etc.).  Because these
various reactions occur quickly, solution equilibria are used to describe the
segments of  the salinity system.  The ion exchange and solution-dissolution
of the slightly soluble salts are computed using standard approaches  such as
the Gibbs phase rule, Debye-Huckel theory, and Gapon equations.  An iterative
successive approximation procedure is used to determine the chemistry within
the soil profile and in the deep percolation.  The procedure involves
 (1) calculation of lime solubility;  (2) calculation of gypsum solubility;

                                     109

-------
(3) undissociated CaS04 and MgSO4 ion pair reaction;  (4) Na+-Ca++, Mg++-Ca++,
and Na+-NH4 ion exchange; and (5) evaluation of CaCO_ dissociation.

     The movement and distribution of water is either read in or produced by
the unsaturated flow program.  For the soluble ions, a mixing cell concept is
used to simulate solute dispersion and movement.  This assumes complete
mixing occurs in each soil-time increment and that molecular diffusion is
negligible.

     Drainage Model (Dynamic Equilibrium).  Flow moving below the crop root
zone is routed to surface receiving waters via this drainage program.  The
model was designed to predict the response of a subsurface drainage system
consisting of parallel, equally spaced tile drains, to deep percolation
inputs.  Principal outputs are the mid-space water table elevation and drain
discharge.  Deep percolation may be developed for the program through card
read input or directly from the unsaturated flow model.  The mathematical
approach is the general Bureau of Reclamation drainage equations.

     The water table shape is approximated as a fourth degree parabola which
uses specific yield to compute instantaneous rises and declines in water
table positions.  Hooghoudt's equivalent depth is used to correct for near-
drain convergence effects not treated by the Dupuit-Forchheimer assumptions.
The Moody solutions are used to approximate these systems.

     Drainage Model (steady-state).  The general modeling system also
includes a steady-state drainage system model based on potential theory to
define stream tubes.  An average deep percolation rate is computed from the
values determined by the unsaturated flow program.  A system of nodes are
defined at which the infinite series mathematical solutions are computed.
Under this approach, the flow velocities and travel times are estimated
thereby yielding an estimate of the time required for a water volume to reach
the drainage system.  These data are determined for use in the saturated
chemistry and drainage effluent chemistry parts of the program.

     Saturated Chemistry Program.  This program predicts the two-dimensional
distribution of saturated flow chemistry below an irrigated field.  The
stream tube flows predicted by the steady-state drainage model along with
chemical soil analysis serve as model input.  Each stream tube is analyzed
separately and is divided into segments of equal volumes.  Flow is accumulated
until it reaches the pore volume of a segment and then is moved by piston
displacement.  After each displacement, the solution is equilibrated with the
solid and exchangeable phases.  Lateral dispersion and diffusion are neglected
since lateral and longitudinal mixing among stream tubes is not considered.
Chemical transformations are identical to those in the unsaturated chemistry
model except denitrification is simulated in the transition zone between
unsaturated and saturated regions.  The denitrification rate is assumed
temperature dependent but not affected by substrate concentration.  Above a
"saturation level" the rates are assumed to be zero-order and below this
level to be first-order.  The aquifer profile is assumed to be homogeneous
until the heterogeneous nature is computed through consideration of
individual stream tubes.
                                     110

-------
     Drain Effluent Prediction.  This program takes the flow rates from the
steady-state drainage analysis and its corresponding chemistry and aggregates
each stream tube into a monthly drainage discharge.  The salt load from the
drainage system is thus determined by mixing and routing the flows in each
stream tube.

GENERAL;

     This modeling system is well documented, verified in field conditions,
and available for broad use.  However, the system itself is far larger than
required by most research groups and far too complex for use by planning
groups.  Many possible uses can be made .using individual segments listed
above.  Computer and data requirements are substantial.  The modeling system
is probably the largest, most exhaustive available as an integrated package.
Code listings are available from the Bureau of Reclamation.
                                      Ill

-------
MODEL:                                F-2

REFERENCES;

Bureau of Reclamation.  1977.  Prediction of mineral quality of irrigation
return flow, Volume III, Simulation model of conjunctive use and water
quality for a river system or basin.  User's Manual.  EPA-600/2-77-179c.
Robert S. Kerr Environmental Research Laboratory, Office of Research and
Development, U.S. Environmental Protection Agency, Ada, Oklahoma.  August.
295 pp.

MODEL SCOPE:

     Simulation is of the water and salinity impacts in a river basin
resulting from urban, industrial, and agricultural water uses.  Formulated as
part of a need for evaluating irrigation return flows in a river system, the
model is primarily oriented to this end.  In addition to the mass balance of
surface and groundwater in a river section, the model simulates reservoir
operations,  and a comparatively complex analysis of the soil salinity system.
Program features both input dumps and error flags.  Intermediate results are
not presented.

INPUT REQUIREMENTS;

     Control parameters are initially read in to completely dictate problem
structures since enough flexibility has been added to allow simulation of
river systems with a multitude of local conditions.  Control parameters
include division of system into segments  (nodes), processes being simulated
in each segment, node structure and labels, and time dimensions.  Input data
descriptive of local hydrochemical processes are input as follows:
(1) aquifer volume and salinity; (2) reservoir conditions and operating
rules; (3)  sediment analyses; (4) irrigated soil chemical and hydraulic
characteristics; (5) hydropower generation parameters; (6) discharge and
chemistry of major surface and subsurface inflows; (7) demands by agricul-
tural, municipal and industrial users; and (8) allocation of flows within an
irrigation system.

SPATIAL AND TIME SCALES;

     Simulation is of hydrochemical systems on a river basin scale or segment
thereof.  Time period utilized is one month at minimum scale to several years
duration of simulation.

COMPUTER CODE STRUCTURE;

     Internal program linkage is accomplished by subroutines and functions
with both common and call statement data transfer (control parameters are
generally passed in the call).  Code is set up to read and write results as
output or onto tape or disk systems.  Approximate core and execution time
requirements are 140 k bytes and 120 central processor seconds, respectively.
Program language is Fortran.


                                     112

-------
BASIC MATHEMATICAL APPROACH;

     The program utilizes an iterative scheme for balancing mass between
"nodes" or subbasins comprising the area being analyzed.  This model includes
substantial simulation of reservoir and stream management.  After initializa-
tion of the system, the various demands are compared with the available
stream flow to determine if a shortage is occurring.  If insufficient water
is available in the river, a series of reservoir checks and manipulations are
made to adjust flow if possible.  Return flow volumes and chemistry are then
computed and checked against initial assumptions.  Iterations are repeated
until mass balance is achieved.

     Of particular interest in this review is the model segments describing
the irrigation return flow system.  Diversions are determined by computing
the consumptive use (using the Penman equation) and then dividing by con-
veyance and farm irrigation efficiencies.  Return flow is then the difference
between demand and consumptive use.  The chemical simulation of the soil
profile is based on a simplified version of the model presented by Dutt, et
al., 1972.

GENERAL:

     Model is very complex in scope and therefore difficult to apply without
some previous experience in hydrosalinity modeling on a large scale.
However, the user manuals are excellent descriptions of the computer system
and can be us^ed to implement the model.  This model has been used extensively
by the Bureau of Reclamation in studies throughout the western United States.
The code is available from the Bureau of Reclamation.
                                     113

-------
MODEL:                                F-3

REFERENCES;

Crawford, N. H. and A. S. Donigian.  1973.  Pesticide transport and runoff
model for agricultural lands.  EPA-660/2-74-013, Office of Research and
Development, U.S. Environmental Protection Agency, Washington, D.C.
December.  317 pp.

Donigian, A.S. and N. H. Crawford.  1976.  Modeling pesticide and nutrients
on agricultural lands.  EPA-600/2-76-043.  Office of Research and Development,
Environmental Research Laboratory, U.S. Environmental Protection Agency,
Athens, Georgia.  February.  211 pp.

MODEL SCOPE:
     Simulation is of runoff, snow accumulations and melt, sediment loss,
pesticide-soil interactions, and soil nutrient transformations.  Analysis
involves transport of pollutants from both surface and subsurface sources.
The model is intended for watersheds up to two square miles in which in-
channel processes and transformations can be neglected.  Parameters simulated
include mass transport of water, sediment, pesticides, nitrogen and phospho-
rus.    The program is set up to operate in either a calibration or production
mode.  In general, the output includes water movements (runoff, infiltration,
evapotranspiration, and storage), sediment transport, pesticide transport and
transformation in various soil layers and in sediment, nutrient transforma-
tions in the soil profile, and leaching of nitrate.  The only salinity
parameter in this version of the model is chloride.  Output can be controlled
to intervals from every 5-15 minutes in the simulation to monthly summaries.

INPUT REQUIREMENTS:

     Data necessary to operate the model can be divided into control,
hydrologic, snow, sediment, pesticides, and nitrogen-phosphorus  parameters.
The control parameters include type of run, units, output interval, snow,
pesticide, nutrient simulation options, input data checks, time intervals for
calculations, watershed structure, and time limits.  For the hydrology sub-
model, input includes soil moisture storage (capacity and initial), length,
slope, roughness, intake, evaporation coefficients, and surface storage of
watershed, groundwater return flow timing, and groundwater distribution and
conditions.  The snow accumulation-melt analysis requires a number of coeffi-
cients describing the energy balance description of the solid-liquid phase
transformations, climatological parameters (radiation, temperature, wind),
and topographic description of the watershed.   To compute the sediment phases
of the model, the user must supply crop cover, tillage effects (depth,
timing, and fine deposits), soil depths and bulk densities, rainfall splash
coefficients, overland flow exponent and wash-off coefficient, and the
initial soil fines deposit.  Pesticides are described using application
timing, solubility, fixing capacity, coefficients in the Fruendlich
adsorption/desorption equation, decay rate, and control coefficients noting
how the pesticide is applied and whether or not a single-valued adsorption/
desorption algorithm is to be used.  Nitrogen and phosphorous simulation


                                     114

-------
begins by first specifying time steps, number of fertilizer applications, and
harvest dates.  Then a series of nitrogen reaction rates for mineralization,
immobilization, nitrification, denitrification, and uptake are input.
Similar coefficients are needed for phosphorus  analysis.  Storage of
nitrogen  (in its various forms) and phosphorus is detailed by describing the
adsorbed, dissolved, absorbed, and gaseous phases.

SPATIAL AND TIME SCALES:
     Calculations are based on a 5 or 15 minute real time interval with
results lumped in hourly, daily or monthly summaries.  The spatial scale is
of a 1-2 square mile agricultural watershed.

COMPUTER CODE STRUCTURE:
     The basic code structure is one of a main program controlling the order
and presentation of input, calculations, and output with subroutines encom-
passing specific submodel tasks.  The code is written in Fortran IV and run
on a IBM 360/67 computer but appears readily adaptable to other machines.
The code is very lengthy and is best operated in compilation and production
steps.  The number of pages of output can be several hundred with even
routine analysis.  Neither core requirements nor execution times are noted,
however, both are comparatively large on visual inspection.

BASIC MATHEMATICAL APPROACH:
     This model is too large for a thorough description of the mathematical
approach.  Consequently, only a brief summary will be presented here with
more detail left to the interested user.  The code involves six primary
subprograms:   (1) MAIN which orders overall execution; (2) LANDS which
simulates the hydrology and snow;  (3) SEDT, the analysis of sediment,
(4) ADSRB, the pesticide adsorption and removal program;   (5) DEGRAD,
pesticide degradation; and  (6) NUTRNT, which is the simulation of nutrient
transformation and removal.

     The LANDS subprogram is derived from the Stanford Watershed Model
presented in various publications.  The hydrologic response of a watershed to
inputs of precipitation and evaporation involves a water budgeting procedure
which accounts for evapotranspiration,.surface storage and retardance,
infiltration, groundwater storage depletion and return flow, and overland
flow.

     The SEDT subprogram is a simulation of sheet and rill erosion involving
the processes of particle detachment by rain drop impact and transport by
overland flow.  The particles detached during an interval are expressed as a
power function of vegetal cover, precipitation, and a detachment coefficient.
Overland transport is expressed as a power function of the detachment
coefficient, initial volume of soil fines, and overland flow rate.  Attempts
have been made to include the effects of tillage operations by allowing the
user to redefine initial volumes of detachable soil fines.  Vegetal cover
effects are specified monthly by the user with linear interpolation during
the month.

                                     115

-------
     The subprogram ADSRB is based on two alternatives.  First, the standard
single-valued Freundlich adsorption/ desorption isotherm is used to separate
adsorbed and dissolved pesticide concentrations.  The second approach
involves a multiple valued adsorption/desorption process taken from the work
of Davidson, et al. (1975) reported elsewhere in this writing.

     The degradation and volatilization of pesticides are simulated in the
subprogram DEGRAD.  Volatilization is ignored in the model's present version
and degradation is assumed to follow a first-order decay function.

     The nutrient simulation model (NUTRNT) attempts to predict nitrogen and
phosphorus  losses due to erosion, surface washoff, leaching, and biological
conversion.  Phosphorus  is assumed to exist as organic phosphorus,  solid
phosphate, dissolved phosphates, and phosphorus  adsorbed by the plants.
Transformations are all based on first-order kinetics  (temperature dependent).
Nitrogen processes (mineralization, immobilization, nitrification, denitri-
fication, and uptake)  are based on first-order kinetics using literature
values of the coefficients.  Ammonium adsorption/desorption is also included
in the model.  Solution of the coupled differential equations is accomplished
with a simple Euler integration scheme.

GENERAL:
     This model is currently being improved and modified for a more general
use.  The nutrient phases remain unverified at this time although the surface
hydrology and pesticide portions are tested against actual field data.  Code
listing is presented in the reference.
                                     116

-------
MODEL:                                F-4

REFERENCES:

Dixon, N.P., D.W. Hendricks, A.L. Huber, and J.M. Bagley.  1970.  Developing
a hydro-quality simulation model.  Report PRWG 67-1.  Utah Water Research
Laboratory, Logan, Utah.  June.  193 pp.

MODEL SCOPE;

     Simulation is of the electrical conductivity, dissolved oxygen,
temperature, BOD, and coliform count in the inflows and outflows of river
subbasins.  Model output includes the spatial profiles of flow rates in a
river section which occur as stream inflow, natural diffuse inflow, ground-
water inflow, surface irrigation return flow, tributary inflow, and municipal-
industrial effluent discharges.  These flows are mixed in a stream segment
and adjusted for evapotranspiration, system storage changes, and precipita-
tion to arrive at the section outflows.  The chemical phases of the model are
based on computing monthly averages except if diurnal distribution of stream
temperature and/or dissolved oxygen are requested.

INPUT REQUIREMENTS;

     The model is divided into four submodels:   (1) hydrology; (2) electrical
conductivity; (3) temperature; and  (4) dissolved oxygen.  Attempts to model
BOD and coliform counts were unsuccessful and not suggested for general
application.  A system control submodel is used to control the order of
calculations, input data, and coordinate submodel interactions.

     Inputs to the hydrology submodel consists of precipitation, measured
stream inflows in the main channels, measured surface inputs, unmeasured
surface and subsurface inflows, cropland diversions, reservoir storage,
municipal-industrial diversions, pumped water, surface exports, and mean air
temperature.  The results from the hydrology submodel detail  the flows
within the river section.  These flows are characterized by input electrical
conductivities or flow versus electrical conductivity relationships.  Simula-
tion of stream temperature requires data describing the temperature of the
respective inflows to a river section along with stream depth, discharge, and
width to depth ratios to calculate heat equilibrium rates.  The dissolved
oxygen submodel requires that the initial oxygen level, BOD^2 concentration,
and rate constants be known.  Flow data are supplied by other submodels.

SPATIAL AND TIME SCALES;

     Although the model has the capacity for day to day change in stream
temperatures and dissolved oxygen, the effective time resolution is one
month.  The scale encompasses major river segments in which uses and natural
inputs tend to be lumped rather than distributed.
                                     117

-------
COMPUTER CODE STRUCTURE;

     The program code is in the form of a single main program divided into
segments by return controls and is coded in Fortran IV for a Univac 1108
computer.  Core requirements listed are 20 k bytes of 32 bit words.  Execution
time is on the order of 30 seconds.

BASIC MATHEMATICAL APPROACH;

     The hydrologic submodel is a mass balance simulation of the water flow
along the stream segment.  Precipitation falling on the land surface is
routed through a snow melt equation when the seasonal transitions occur.
Consumptive use is estimated with the Soil Conservation Service version of
the Blaney-Criddle equation.  The transport of water into and out of the
groundwater is based on simple delay coefficients, and average storage
values.  This groundwater segment of the hydrologic submodel is actually
evaluated by residual calculations by adjusting the delay coefficients until
reasonable flows are predicted.

     The basis of the salinity submodel is derived from functions relating
discharge to salinity.  Such expressions are developed for each flow segment
used as input and then computed for the output values based on flow mixing.

     Temperatures of the surface flows are predicted using linear regressions
of air temperature and water temperature. Respective flows given as input are
mixed to predict outflow temperatures.  Some simulation is possible using
sine functions of air temperature.  The changes in water temperature inside a
stream segment is based on a LaBouquet decay function.

     The dissolved oxygen submodel utilizes the Streeter-Phelps oxygen sag
equation.  The BOD impact on dissolved oxygen is evaluated by the Hansen-
Frankel equation.  The first-order rate coefficients are taken from the
literature.

GENERAL;

     The reference cited above contains both a listing and a detail set of
user instructions.  The hydrology submodel has been substantially updated in
other work at Utah State University.
                                     118

-------
MODEL;                                F-5

REFERENCES:

Hill, R.W., E.K. Israelsen, and J.P. Riley, 1973.  Computer simulation of the
hydrology and salinity system within the Bear River Basin.  PRWG 104-1.
Utah Water Research Laboratory, Utah State University, Logan, Utah.  122 pp.
Narasimhan, V.A., and E.K. Israelsen, 1975.  A water-land use management
model for the Sevier River Basin, Phases I and II.  Information Bulletin 24.
Utah Water Research Laboratory, Utah State University, Logan, Utah.
September.  44 pp.
Huber, A.L., E.K. Israelsen, R.W. Hill, and J.P. Riley, 1976.  Basin
simulation assessment model documentation and users manual.  PRWG 201-1.
Utah Water Research Laboratory, Utah State University, Logan, Utah.  August.
30 pp.

MODEL SCOPE:

     Simulation involves the basic hydrologic relationships in river basins
and/or watersheds.  The model incorporates relationships describing hydrologic
processes which are linked together by the conservation of mass principle.
Salinity is added to simulation by assuming the hydrologic processes change
the chemical concentration in the system by storing, concentrating, diluting
and/or picking up salts.  Salinity related processes are specifically
developed for irrigated soils and small storage reservoirs.
                                  f1
     The model concept is one of inflow-outflow = change in storage.  Inflows
simulated include stream inflows, tributary inflow, precipitation, subsurface
inflow, and imports.  Outputs are surface and subsurface flows, evaporation
and evapotranspiration, and exports.  Internal system- simulations consider
canal and groundwater diversions, surface and subsurface irrigation return
flow, and soil moisture groundwater systems.  Reservoir operations are
provided with a specific submodel.

     Use of this model falls under two categories.  First, given the
necessary input data, the simulation of various hydro-salinity flows is
compared with actual measured values in order to calibrate the model.  A
pattern-search optimization technique is utilized to systematically optimize
the calibration.  The second model use which is based on proper calibration
is to test management alternatives  (changes in cropping patterns, alternative
reservoir operations, and improved irrigation practices) on the water quality
of the system outflows.

INPUT DATA:

     Control parameters allow user to select one of four input options:
(1) read run data for reduction and handling internally;  (2)  read data
directly to simulation system;  (3) read data from computer memory; and
(4) read parameter bounds for calibration process.  Other control parameters
set boundaries for calibration process.  Input of a hydrologic nature which
must be supplied includes monthly percentage of daylight hours, crop and

                                     119

-------
phreatophyte consumptive use coefficients, canal diversion, land use acreages,
precipitation, soil moisture characteristics, temperature, river inflow,
tributary inflow, subsurface inflow, surface outflow, and reservoir operations.
Salinity data are also necessary for the above noted water flows.

SPATIAL AND TIME SCALES:

     Simulation is of single or multiple river system subbasins.  Analysis is
partially dynamic in nature although steady-state is assumed within a time
frame of 1 month.  Time resolutions are monthly with multiple year analyses
inherent.  Model is two-dimensional in nature.

COMPUTER CODE STRUCTURE:

     Internal program linkage is accomplished by subroutines with call and
common data transfer.  The program in its initial versions was written and
utilized on a digital-electric analog hybrid computer at Utah State University.
Hydrology is currently available in a digital format.  The digital segments
should be adaptable to other computer systems without much difficulty.
Language is Fortran IV.

     Execution time and memory requirements are not given.  General
examination indicates core requirements on the order of 50-70 k bytes and
execution times on the order of 3-5 minutes for major calibration-simulation
runs.

BASIC MATHEMATICAL APPROACH;

     The basic mathematical relationship in the model is the conservation of
mass within a river basin subarea:

          (PPT + QSI + QGI) -  (GET + QSO + QGO) = AS

where,    PPT = precipitation
          QSI = total surface inflow
          QGI = total groundwater inflow
          GET = evapotranspiration
          QSO = total surface outflow
          QGO = total groundwater outflow
           AS = change in system storage in snow, soil moisture, surface
                reservoirs, and groundwater.

     The principal model variable is QSO.  As a generally known quantity, the
model simulations of QSO are compared with measured values.  Calibration
involves adjusting various time lag and simulation coefficients until the
error in the predicted values of QSO is minimal.  The minimization approach
has been termed a pattern-search technique which is essentially a process of
iteratively varying the parameters to achieve minimum variance.
                                      120

-------
GENERAL;

     This model is one of the most recent in a long-term program of
hydrologic modeling at Utah State University.  Many of the earlier models or
versions thereof have been incorporated or omitted from this model based on
applications in a number of field situations.  A listing is included in the
reference by Huber, et al. (1976).
                                     121

-------
MODEL;                                F-6

REFERENCES:

Walker, W.R.  1970.  Hydrosalinity model of the Grand Valley.  Unpublished
M.S. Thesis, Civil Engineering Department, Colorado State University, Fort
Collins, Colorado.  August.  94 pp.

MODEL SCOPE;

     This simulation is limited to the hydrologic and salinity systems in an
irrigated valley.  Diversions from reservoirs, rivers, tributary inflows, or
groundwater are routed through the agricultural system.  Surface as well as
subsurface flows and flow qualities are predicted along with consumptive use.

INPUT REQUIREMENTS;

     Input data consist of mean monthly values for temperature, percent
daylight hours, Blaney-Criddle climatic coefficients, areas of each land use,
Blaney-Criddle crop growth stage coefficients, river inflows, tributary
inflows, river outflows, lateral diversions, surface drainage outflows,
drainage base flows, soil moisture storage capacities, root zone diversions,
exports, imports, water table fluctuations, precipitation, and equilibrium
salinity concentrations of each of these flows.  A series of groundwater data
are also required.  These include the number of strata at the groundwater
basin's return flow point, and the number of normal locations where the
gradient in each of these strata is defined.  At the outflow point or close
by, an initial value of hydraulic conductivity and cross-sectional area is
required.

SPATIAL AND TIME SCALES;

     The scale of the model is a macroscopic simulation of an irrigated area
and the hydrology immediately surrounding it.  A time resolution of one month
is taken.  Analysis is steady-state.

COMPUTER CODE STRUCTURE;

     The computer code consists of a main program and subroutine served by
both common and call statement data transfer.  Written in Fortran IV, it
should be adaptable in most computers.  Core requirements are approximately
43 k bytes and time needs are about 25-30 central processor seconds per run.

BASIC MATHEMATICAL APPROACH;

     A mass balance is made of the water flow system and the salt system is
attached by multiplying by the concentrations associated with each segment of
the system.  Flows entering the main conveyance are segregated into seepage
by a conveyance efficiency, lateral diversions as input, and operational
wastes by difference.  Lateral diversions are likewise delineated into
seepage, root zone diversions, and tailwater.  Finally, root zone diversions
are allocated to root zone storage, deep percolation or consumptive use.

                                     122

-------
Summing the surface and subsurface flows from each breakdown yields an
estimate of the respective return flows.  At this point an estimate of the
groundwater return flows is made on the basis of a groundwater model.   A
comparison of the mass balance and groundwater model estimates of subsurface
return flows is made in a manner which calculates an adjusted aquifer
hydraulic conductivity by month.  The model is then systematically adjusted
until the values of hydraulic conductivity are homogeneous.  When the water
system is satisfactorily modeled, the salinity system is added.

GENERAL:

     The model is a second or third generation modeling concept initiated
at Utah State University.  Consequently, these early models will not be
described; however, this model is an approach that can be applied with
generally available data.  The evapotranspiration and root zone analysis
needs updating.  Application has been made with good accuracy in the
Grand Valley of western Colorado.
                                     123

-------
                                   TECHNICAL REPORT DATA
                            (Please tsad launtctions on the reverse before completing)
1. REPORT NO.
  EPA-600/2-78-144
                              2.
                                                           3. RECIPIENT'S ACCESSIOWNO.
4. TITLE AND SUBTITLE

 IDENTIFICATION AND INITIAL  EVALUATION OF
 IRRIGATION RETURN FLOW MODELS
                                   5. REPORT DATE
                                    July 1978 issuing  date
                                   6. PERFORMING ORGANIZATION CODE
7. AUTHORIS)

 Wynn R.  Walker
                                                           8. PERFORMING ORGANIZATION REPORT NO.
               \NIZATION NAME AND ADDRESS
9. PERFORMING ORGANIZATION NAME AND
 Irrigation Hydrology Company
 P.  0.  Box 1544
 Fort Collins, Colorado  80522
                                   10. PROGRAM ELEMENT NO.

                                    1HB617
                                   11. CONTRACT/GRANT NO.

                                    R-804515
 12. SPONSORING AGENCY NAME AND ADDRESS
 Robert S.  Kerr Environmental  Research Laboratory-Ada,OK
 Office of  Research and Development
 U.S.  Environmental Protection Agency
 Ada,  OK  74820
                                    13. TYPE OF REPORT AND PERIOD COVERED
                                    Final
                                   14. SPONSORING AGENCY CODE
                                    EPA/600/15
 15. SUPPLEMENTARY NOTES
 A oroad based literature review was undertaken to  identify studies that had yielded
 digital computer models applicable to irrigation return flow (IRF) systems.   The
 programs not listed in technical reports or papers were requested from the various
 authors.   The results of this  work are 43 computer models applicable all or in part
 to the  analysis of IRF's and their quality.  A brief evaluation of each model is given

 IRF modeling technology is well developed theoretically but not completely verified
 due to  the large scale of the  irrigation system.   Most models remain in the research
 sphere  and need to be redefined for the wider utilization of planners.  Field data
 are generally not available to satisfy the input requirements of most IRF  models.
 Accuracies of predictions need to be determined against standardized conditions in
 order to further model development and parameter sensitivities should be investigated
 to isolate the most important  field data.
17.
                                KEY WORDS AND DOCUMENT ANALYSIS
                  DESCRIPTORS
                                              b.IDENTIFIERS/OPEN ENDED TERMS
                                                 c. COSATI Field/Group
 Mathematical Models
 Irrigation
 Water Resources
 Water Pollution
 Water
 Simulation
 Prediction
Water Quality
Soil Water
Soil Chemistry
Evapotranspiration
Irrigation Return Flow
                           72/E
                           68/D
18. DISTRIBUTION STATEMENT

 RELEASE TO PUBLIC
                      19. SECURITY CLASS (ThisReport)
                       Unclassified
                          21. NO. OF PAGES
                           130
                      20. SECURITY CLASS (This page)
                       Unclassified
                                                 22. PRICE
EPA Form 2220-1 (9-73)
                                            124
                                                                     *U.S. BMHMBITPWraieOFFICE BJS— 757-140/1404

-------