
W END OF A A1
Cn
LAKE APOPKA
Figure 70. Lateral view of a portion of crosssection AA» (Figure
23)

c/i
oo
I 2
HEAD INCREMENT, ft
Figure 71. Predicted recharge velocities resulting from various increases
in water table elevation.

Additional information was sought by extending the simulation period.
The predicted contours following 68.8 days of input are given in Figure
72. In the right hand portion of the figure, the advancing high con
centration (10"*) front has been slowed down. This is due to the low
permeability section in the vicinity (see Figure 45).
The high vertical velocities in the upper nodes of the sand formations
dictated that a time step no greater than one day be used. For this
reason, simulation over long time periods would have been very costly.
SUMMARY
The results obtained indicate that it is feasible to model a conjunc
tive surfacegroundwater system. Current data limitations precluded
an adequate calibration of the model. However, the formulation can
be used as a guide for data collection.
While it is desirable to conduct modeling from a systems viewpoint,
it appears that a close integration of the various individual pro
grams would be premature. In their present relatively independent
status, the programs are more accessible for the introduction of im
provements. Also, the high computation cost and programming complex
ities inherent in a fully interfaced approach seem to be unwarranted
in the light of data deficiencies.
SUGGESTIONS FOR FURTHER WORK
Discussed below are possible improvements that might be incorporated
into the methodology and some extensions of its usage. Considered
first are improvements of the individual models. Data collection and
parameter optimization factors are discussed next. Finally, several
comments about the study area are included. It is likely that some
of these comments will be applicable to other locations.
In the formulation of improvements, a compromise must be established
between accuracy and cost. The desirable position of this balance
depends strongly on the nature of the study area and the planned use
of the model. The numerous costefficiency tradeoffs can perhaps be
most effectively evaluated by decisionmakers, who are generally the
ultimate users of models (Kisiel, 1971).
Individual Models
Surface Water Models
Kinematic wave model The little apparent significance of overland
flow in the Lake Apopka basin and the lack of appropriate data precluded
adequate verification of the kinematic wave relationships. Nevertheless,
it is felt that several improvements therein should be given considera
tion, as follows:
159

Q.
LJ
O
O
CO
<
LJ
(T
O
a:
LU
00
O
a:
2.
5.
6.
7 
in •
1C
lo
14
1C .
IO
\l
ill!
ll\V
^
sw
s
J/s
'^
SAN
1
DY S
X
50ILJ
3 ~\
/
•^^K — ••
— — —
y
/
JVKE
AP(
)PKA
X
Z
^
<
Ml
^.
^
\vv
N
R"*"
IE S
FARI
^T"
^;'.
ANm
*s
y?
^/
r so
=^
\
\
\
\
LS/
I
\
S
\\)
All
i
^
•
^
//
/i
?'
Symbols Concentrations, mg/l
IO'1
ID'2
IO'3
ID'4
T^
2
3
4
6
7
8
9
12
14
15
16
• 17
• 18
19
?r
89 10 II 12 13 14 15 16
COLUMN NUMBER
17 18 19
Figure 72. Predicted nitrate concentration contours in crosssection AA!
(Figure 23) after 68.8 days (1656 hours) of input.

1. Allowance should be made for a given catchment to provide flows
to another catchment. Recent reports from the Massachusetts institute
of Technology (Harley, Perkins and Eagleson, 1970) contain useful
information on this matter.
2. A more rigorous interfacing between overland flow and flow through
the unsaturated zone might be desirable. Smith and Woolhiser (1971)
developed a numerical formulation for the simultaneous consideration
of these two regimes. Their methodology, Which is quite general,
includes provision for cascading catchment planes. However, its
computational costs might be relatively high.
3. Attempts should be made to account for the effects of agricultural
practices in the model. As indicated by Hills (1971), the infiltra
tion capacity of a soil depends more on treatment factors (e.g., com
paction) than on its natural characteristics. In 1971, Seginer devel
oped a model for surface drainage of cultivated fields, which gives
consideration to furrows and their hydraulic characteristics. The
formulation developed by the Agricultural Research Service (1971) is
also quite valuable. Testing of Manning coefficient relationships for
vegetated areas (Corps of Engineers, 1956) would be desirable.
4. Much remains to be accomplished in the matter of suspended solids
and phosphorous fluxes in overland flow. A better representation of
erosion phenomena might be obtained by utilizing a Lagrangian modeling
approach, in which the path of the particles are followed as they move
through the system. Additional agricultural data are required to
improve estimates of phosphorous content in suspended soil particles.
Stormwater Management Model The stormwater model offers advantages
in its inherent structure for designating catchments and reading in
data in an organized fashion. Its structure is also well suited for
interfacing with the receiving model. However, several modifications
to the formulation should be contemplated:
1. At present, the model contains calculations for hydraulic elements
which are not found in agricultural areas, such as pipes and gutters.
These calculations were bypassed in this study, but for more intensive
use of the model it might be worthwhile to carry out a more thorough
elimination and rewriting process. The same argument applies to water
quality computations, s*uch as those for biochemical oxygen demand
(BOD) and coliform bacteria counts.
2. Consideration should be given to improvements regarding closer
interfacing with unsaturated zone flows, accounting for agricultural
practices and sediment load calculations. The same specific comments
provided for the kinematic wave model apply herein.
161

Receiving water model— The performance of the flow portion of this
model, in both the drainage canals and Lake Apopka, was satisfactory.
Due to the flexibility of its linknode structure and its ability to
consider various boundary conditions, the model has great applica
bility to practical situations. Its major drawback lies in its high
computational cost. Suggested improvements lie in the area of quality
calculations:
1. For a complete understanding of nutrient behavior, it would be
desirable to incorporate an ecosystem model in the quality subroutines.
Such models have been proposed by numerous investigators, including
Chen (1969), Bloom, Levin and Raines (1969) and Chen and Orlob (1972).
Since the amount of data required to calibrate an ecosystem model is
rather massive, the development of a data collection strategy would
appear worthwhile. In this regard, Park and Wilkerson (1971) provide
insight into the interaction between statistical studies and sampling
programs. In a recent paper, Moore (1972) discussed the use of filter
theory in data collection, and the tradeoffs between spatial and tem
poral frequency in sampling.
2. Since often appreciable amounts of nutrients are contained in soil
particles, it would be pertinent to devote further attention to the
phenomenon of sediment transport. The consideration of this phenome
non within the receiving bodies of the Lake Apopka basin (and other
agricultural basins) would appear to be of interest. In a study of
cohesive soils, Partheniades (1965) concluded that erosion rates were
rather independent of the concentration of suspended sediment, and
arrived at several other generalizations regarding erosion and deposi
tion factors. In 1968, Etter <5t aJL discussed the importance of parti
cle sizes and their relationship to flocculation. In addition to in
formation contained in such studies, data are required on nutrient up
take and realease rates from sediments.
Unsaturated Zone Model—
The most important suggested improvements for this model are on the
subject of water quality. Following is a list of suggestions:
1. It might be desirable to integrate this model more fully with an
overland flow formulation (Smith and Woolhiser, 1971), as discussed
previously.
2. There is a general scarcity of data on hydraulic conductivities
under unsaturated conditions, especially for muck soils. \ More infor
mation on this matter would be desirable.
3. Much attention should be paid to the impact of agricultural activi
ties on flow conditions in soils. Some of the basic relationships
developed in this study for water uptake by plants could be improved in
162

their representation of seasonal patterns and the effects of different
crops. Another possible feature to be included would be the flow of
water through the voids between the plant roots and the soil. This and
other features are contained in a model developed recently by the
Agricultural Research Service (ARS) (1971)•
4. A considerable amount of work remains to be done on the determina
tion of source/sink rates for nutrients in soils. At present, the
qualitative understanding of the mechanisms affecting these rates is
very limited. It is suggested that the model developed herein be used
as a guide for data collection and for verification purposes. It is
also thought that the ARS model could be improved by the incorporation
of water quality routines.
Groundwater Model—
It is felt that the decision to utilize a steadystate twodimensional
model, explained in previous sections, was justified in terms of cost
and the accuracy desired for most agricultural situations. However,
in basins where pumpage and other boundary conditions produce signifi
cant changes in piezometric conditions over time and space, the use of
an unsteadystate threedimensional formulation (Freeze, 1971) may be
necessary. Most of the remaining comments deal with general approaches
to model calibration:
1. In groundwater systems, it is commonplace to find that more than
one set of heads and permeabilities will satisfy certain known flow
conditions. This problem of indeterminacy poses serious difficul
ties to the calibration of both the flow and quality portions of the1
model. The possible use of an optimization approach to the estimation
of model parameters is discussed later.
2. It appears worthwhile, on a longrange basis, to pursue the global
energy approach proposed by Back and Hanshaw (1970). The success of
this method, which was discussed previously, would depend on both a
better understanding of energy phenomena and intensive sampling. It
would seem appropriate to devote the original modeling efforts in
connection with this approach to welldefined geologic systems and
substances of simple chemical behavior.
3. Additional quantitative information should be obtained on source/
sink reactions of nutrients in groundwater. Since agricultural
phosphorous is likely to be fixed to the shallow soil zone, more
attention might be devoted to anaerobic nitrogen reactions.
Data Collection and Parameter Estimation
The problem of indeterminacy, which is quite significant in ground
water models, is also present in other subsystem formulations. Thus
163

the consideration of data collection aspects is extremely complex,
Several data collection strategies have been developed whose basic
objective is to minimize the variance between observed and predicted
values (Matalas, 1967). However, in the modeling of a conjunctive
surfacegroundwater system the precise definition of an objective func
tion is an elusive matter. This problem is due to the stepwise (flow
followed by quality) nature of the predictive process, and the large
number of subsystem models. In designing sampling programs, it would
seem reasonable to assign higher priority to the data that are more
likely to affect pollution abatement decisions. Kisiel (1971) has
provided useful insight into the use of decision theory to determine
the value of data.
In view of the imperfect nature of any data base, it is also worth
while to enhance the predictive quality of models by improving the
estimates of various parameters. Optimization methods for flow resis
tance in overland flow (Schreiber and Bender, 1972) and for per
meability in groundwater (Kleinecke, 1971) are available. The latter
contains restrictive assumptions such as homogeneity of the medium
and the occurrence of flow in two dimensions only. Since available
optimization techniques involve only flow phenomena, the development
of counterparts would be required to calibrate quality processes,
The multiplicity of subsystems within any given basin suggests the use
of a decomposition, rather than integrated, approach for both data col
lection and parameter optimization. Much remains to be accomplished
on the quantitative understanding of interactions between subsystem
models through boundary conditions, and on the effect of flow predic
tion improvements upon water quality modeling. Regarding the latter,
perhaps the global energy approach suggested by Back and Hanshaw (1970)
could be used to eliminate the need for the classical twostep (flow,
then quality) approach.
Suggestions for the Lake Apopka Basin
Due to the current interest in the water quality of Lake Apopka and
its possible restoration, it seems pertinent to summarize the main
data needs and potential model uses:
1. Available data were inadequate to determine a reasonably accurate
hydrologic budget of the muck farms area (Heaney, Perez and Fox', 1972).
Pumpage figures are quite rough and no records are kept on the opera
tion of the flap valves. It is suggested that better data be obtained
for these flows. Groundwater flows in the area also merit further
attention.
2. In the muck farms, the impact of fertilizer application appears to
be less significant than that of nutrient releases by the soil. Data
164

leading to a better understanding of the mechanisms affecting this
release would be very desirable.
3. The nitrogen contributions from the citrus groves to the groundwater
system might be appreciable. However, due to the piezometric conditions,
any fluxes into Lake Apopka are likely to be small. The only signi
ficant contribution to the lake could occur at Gourd Neck Spring.
4. The results of the receiving model indicate that the horizontal
movement of substances in Lake Apopka is quite slow. Thus, it might
be worthwhile to incorporate the phenomenon of diffusion into the
model. A better quantitative understanding of source/sink reactions
between the lake water and the bottom sediments would be very help
ful. Research on sediment transport aspects, including particle size
and scouring characteristics, should be useful. Also, the mixing effect
of wind should be included. This effect might account for the fairly
uniform (in the horizontal sense) pattern of nutrient concentrations
in the lake.
5. Overland flow appears to be of very limited significance in the
Lake Apopka basin. In future modeling efforts in this basin, it
might be desirable to bypass consideration of this phenomenon, except
in the case of very high intensity storms,
6. The methodology developed herein could be used as a guide for
data collection in this and other basins. Following satisfactory
verification, the model could also be utilized as a tool for simu
lating different pollution abatement measures.
165

SECTION X
REFERENCES
1. Agricultural Research Service, U. S. Dept, of Agriculture, Techni
cal Bulletin No. 1435, Beltsville, Maryland. 1971. (unpublished).
2. Anderson, W. and B. F. Joyner. Availability and Quality of Surface
Water in Orange County, Florida. Map Series No. 24, Florida Board
of Conservation, Tallahassee, Florida. 1969.
3. Anderson, W. Hydrologic Considerations in Draining Lake Apopka—
A Preliminary Analysis. Open File Report 71004, U. S. Geological
Survey, Tallahassee, Florida. 1971.
4. Ashcroft, G., et_ al_. Numerical Method for Solving the Diffusion
Equation: I. Horizontal Flow in SemiInfinite Media. Soil Sci
Am Proc. 2(5:522525, 1962.
5. Bachmat, Y. and J. Bear. The General Equations of Hydrodynamic
Dispersion in Homogeneous, Isotropic Porous Mediums. Journal
of Geophysical Research. £9(12):25612567, June 1964.
6. Back, W. and B. B. Hanshaw. Comparison of Chemical Hydrogeology
of the Carbonate Peninsulas of Florida and Yucatan. Journal of
Hydrology. j_0:330368, 1970.
7. Bailey, G. W. Role of soils and Sediment in Water Pollution
Control. EPA Southeast Water Lab, Athens, Georgia, March 1968.
8. Bates, T. E. and S. L. Tisdale. The Movement of NitrateNitrogen
through Columns of CoarseTextured Material. Soil Sci Soc Am Proc.
Vol.21, 1957. ' " :
9. Bear, J. On the Tensor Form of Dispersion in Porous Media.
Journal of Geophysical Research. 60_(4) : 11851198, Apr.il 1961.
10. Bear, J. Dynamics of Fluids in Porous Media. AmericanElsevier.
1972.
166

11. Bedient, P. B. A TwoDimensional Transient Numerical Model for
Radionucleide Transport in Tidal Waters. Master of Science Thesis,
Environmental Engineering Dept., University of Florida, Gaines
ville. March 1972.
12. Benson, N. and R. N. Barnette. Leaching Studies with Various
Sources of Nitrogen. J Amer Soc Agron. Vol. 31, 1939
13. Bigger, J W., D. R. Nielson and K. K. Tanji. Comparison of
Computed and Experimentally Measured Ion Concentrations in Soil
Column Effluents. Amer Soc Agric Engr Trans. £: 784787.
14. Blaney, H. F. and W. D. Griddle. Determining Water Requirements
in Irrigated Areas from Cliraatological and Irrigation Data.
U. S. Dept. of Agric. Technical Paper No. 96. 1950.
15. Bloom, S. G., A. A. Levin and G. E. Raines. Mathematical Simula
tion of Ecosystems A Preliminary Model Applied to a Lotic Fresh
water Environment. Battelle Memorial Institute Report, Columbus,
Ohio. April 1969.
16. Bravo, C. A. et. al^. A Linear Distributed Model of Catchment
Runoff. MIT Hydrodynamics Lab Report No. 125, Cambridge, Mass.
June 1970.
17. Bredehoeft, J. D. and G. F. Finder. Digital Analysis of Areal
Flow in Multiaquifer Groundwater Systems: A Quasi ThreeDimen
sional Model. Water Resources Research. £(3):883888, 1970.
18. Brezonik, P. L. Eutrophication: The Process and its Modeling
Potential. In: Modeling the Eutrophication Process. Proc.
of St. Petersburg Workshop, Department of Environmental Engineer
ing, University of Florida, Gainesville. 1969.
19. Brutsaert, W. The Permeability of a Porous Medium from Certain
Probability Laws for Pore Size Distribution. Water Resources
Research, £:425434, April 1968.
20. Butson, D. K. and G. M. Prine. Weekly Rainfall Frequencies in
Florida. Circular S187. Institute of Food and Agricultural
Science (IFAS), University of Florida, Gainesville. April 1968.
21. Champlin, J. B. F. and G. G. Eichholz. The Movement of Radio
active Sodium and Ruthenium through a Simulated Aquifer. Water
Resources Research. 4(1):147158, Feb. 1968.
167

22. Chen, C. W. Concepts and Utilities of an Ecologic Model. In:
Modeling the Eutrophication Process. Proc. of St. Petersburg
Workshop, Department of Enviromental Engineering, University of
Florida, Gainesville. November 1969. p. 124126.
23. Chen, C. W. and G. T. Orlob. Ecologic Simulation for Aquatic
Enviroments. Water Resources Eng., Inc., Walnut Creek, Calif.
N.T.I.S. No. PB218 828, Dec. 1972. 156 pp.
24. Choate, R. E. and D. S. Harrison. Irrigation by the Accounting
Method. Agric Engr Mime.o Report No. 685, University of Florida,
Gainesville. May 1968.
25. Crawford, N, H. and R. K. Linsley. Digital Simulation in Hydro
logy: Stanford Watershed Model IV. Technical Report No. 39,
Stanford University. July 1966.
26. Dagan, G. Perturbation Solutions of the Dispersion Equation in
Porous Mediums. Water Resources Research. £(1):135142, Feb.
1971.
27, Day, P. R. and W. M. Forsythe. Hydrodynamic Dispersion of Solutes
in the Soil Moisture Stream. Soil Sci Soc Amer Proc. 21;477
480, 1957.
28. Diskin, M. H. On the Computer Evaluation of Thiessen Weights.
Journal of Hydrology. 11(1):6978, July 1970.
29. Domenico, P. A., D. A. Stephenson and G. B. Maxey. Groundwater
in Las Vegas Valley. Technical Report No. 7, Desert Research
Institute, Reno, Nevada. April 1964,
30. Domenico, P. A. Concepts and Models in Groundwater Hydrology.
New York, McGrawHill, 1972.
31, Dyer, K. L. Unsaturated Flow Phenomena in .Panoche Sandy Clay
Loam as Indicated by Leaching of Chloride and Nitrate Ions.
Soil Sci Soc Am Proc. 29^121126, 1965.
32. Dyer, K.L. Interpretations of Chloride and Nitrate Ion Distri
bution Patterns in Adjacent Irrigated and NonIrrigated Panoche
Soils. Soil Sci Soc Amer Proc. 29:170176, 1965.
33. Eagleson, P. S. and W. 0. Maddaus. A Distributed Linear Representa
tion of Surface Runoff. MIT Hydro. Lab Report No. 115. June 1969.
34. Eagleson, P. S. Dynamic Hydrology. New York, McGrawHill, 1970.
168

35, Edinger> J. E. and J. E. Geyer. Heat Exchange in the Environ
ment. Edison Electric Institute, Publication 65902. June 1965.
36. Einstein, H. A. Part II. River Sedimentation, In: Handbook of
Applied Hydrology, Chow, V. T. (ed.). New York, McGrawHill,
1964. p. 1752.
37. EngineeringScience, Inc. and Grunwald, Crawford and Associates.
Central Fresno County Water and Liquid Wastes. Vol. 1, March
1970. p. VI1 to VI8.
38. Environmental Protection Agency. Storm Water Management Model.
(4 volumes), EPA, Washington, D.C., Reports 11024 DOC07/71,
08/71, 09/71, 10/71, September 1971.
39. Etter, R. J., R. P. Hoyer, E. Partheniades, and J. F. Kennedy.
Despositional Behavior of Kaolinite in Turbulent Flow. Journal
of the Hydraulics Division, ASCE. 9£(HY6):14391452, Nov. 1968.
40. Federick, J. C. Solving Disposal Problems in Unsewered Areas.
Sewage Work Eng. Vol. 19, 1948.
41. Florida State Board of Health. Summary Report of Lake Apopka.
Jacksonville, Florida. 1964.
42. Fok, Y. S. and V. E. Hansen. OneDimensional Infiltration into
Homogeneous Soil. J Irrig and Drainage Div, ASCE. £2_: 3546,
Sept. 1966.
43. Fok, Y. S. OneDimensional Infiltration into Layered Soils.
J Irrig and Drainage Div, ASCE. 96_: 121129, June 1970.
44. Forbes, R. B. Water Quality Studies, Zellwood Drainage and Water
Control District. Proc Soil and Crop Science Soc of Florida.
2£:4248, 1968.
45. Forbes, R. B. Letter to P. B. Rhoads of the East Central Florida
Regional Planning Council, dated Sept. 13, 1971.
46. Fox, A. and D. J. Allee. Economic Evaluation of External Effects
of Fertilizer Use. Inr Relationship of Agriculture to Soil and
Water Pollution. Proc Cornell University Conf on Agricultural
Waste Management, Rochester, New York. Jan 1970. p. 185191.
47. Free, G. R. E., G. M. Browning and G. W. Musgrave. Relative Infil
trations and Related Physical Characteristics of Certain Soils.
U. S. Dept. of Agric. Bulletin No. 729, 1940.
169

48. Freeze, R. A. and P. A. Witherspoon. Theoretical Analysis of
Regional Groundwater Flow. 1. Analytical and Numerical Solutions
to the Mathematical Model. Water Resources Research. 2:641
656, 1966.
49. Freeze, R. A. and P. A. Witherspoon. Theoretical Analysis of
Regional Groundwater Flow 2, Effect of Water Table Configuration
and Subsurface Permeability Variation. Wat er jtes ources Res earch.
3^:623624, 1967.
50. Freeze, R. A. Program POTEV  Potential Evapotranspiration
Calculations by the Thornthwaite and Penman Methods; Holmes
and Robertson Moisture Budget Technique. Inland Waters Branch,
Dept. of Energy Mines and Resources, Calgary, Alberta, Canada.
Jan. 1967.
51. Freeze, R. A. The Mechanism of Natural Groundwater Recharge and
Discharge. 1. OneDimensional, Vertical, Unsteady, Unsaturated
Flow above a Recharging or Discharging Groundwater Flow System.
Water Resources Research. £(1):153170, Feb. 1969.
52. Freeze, R. A. and J. Banner. The Mechanism of Natural Ground
Water Recharge and Discharge. 2. Laboratory Column Experiments
and Field Measurements. Water Resources Research. £(1) :138155
Feb. 1970.
53. Freeze, R. A. ThreeDimensional, Transient, SaturatedUnsaturated
Flow in a Groundwater Basin. Water Resources Research. ^(2):
347366, 1971.
54. Gambell, A. W. and C. W, Fisher. Occurrence of Nitrate and Sul
fate in Rainfall. Journal of Geophysical Research. Vol. 69, 1964.
55. Gburek, W. J. and W. R. Heald. Effects of Direct Runoff from
Agricultural Land on the Water Quality of Small Streams. In:
Relationship of Agriculture to Soil and Water Pollution. Proc
Cornell University Conf on Agricultural Waste Management, Roches
ter, New York. Jan. 1970. p. 6168.
56. Gilchrist, A. N. and A. G. Gillingham. Phosphate Movement in
Surface Runoff Water. New Zealand Journal of Agric Research.
.13:225231, 1970.
57. Gottschalk, L. C. Part I. Reservoir Sedimentation, In: Hand
book of Applied Hydrology, Chow, W. T. (ed.)« New York, McGraw
Hill, 1964, p. 1724.
170

58. Glymph, L. M., H. N. Holtan, and C. B. England. Hydrologic
Response of Watershed to Land Use Management. Jour ofIrrig
and Drainage Division, ASCE. 97_(IR2):305318, June 1971.
59. Green, D. W., et_ al^. Numerical Modeling of Unsaturated Ground
water Flow and Comparison of the Model to a Field Experiment.
Water Resources Research. 6^(3):862874, June 1970.
60. Hadley, G. Linear Algebra. Reading, Massachusetts, Addison
Wesley Co., 1964.
61. Haimes, Y. Y., R. L. Perrine, and D. A. Wismer. Identification of
Aquifer Parameters by Decomposition and Multilevel Optimization.
Water Resources Center Publ. No. 123, University of California,
Los Angeles, March, 1968.
62. Hanks, R. J. and S. A. Bowers. Numerical Solution of the Mois
ture Flow Equation for Infiltration in Layered Soils. Soil Sci
Soc Am Proc. J26.: 530534, 1962.
63. Hantush, M. S. and C. E. Jacob. NonSteady Flow in an Infinite
Leaky Aquifer. Trans Amer Geophysical Union. 36^:95100, 1955.
64. Harley, B. M., F. E. Perkins, and P. S. Eagleson. A Modular
Distributed Model of Catchment Dynamics. MIT Hydrodynamics
Lab Report No. 133, Cambridge, Massachusetts. Dec. 1970
65. Harrison, D. and H. A. Weaver. Some Drainage Characteristics
of a Cultivated Soil in the Everglades. Soil and Crop Science
Soc of Florida. 18}184192, 1958.
66. Harrison, D. S. Applying Basic Soil Water Data to Water Control
Programs in Everglades Peaty Muck. ARS 4140, Agric. Res.
Service, USDA, No. 1960.
67. Harrison, D. S., E. H. Stewart and W. H. Spier. Some Considera
tions for Drainage of Flatwoods Soils Used in Vegetable Produc
tion. Proc Florida State Horticultural Society. Vol. 73, 1960.
68. Hashemi, F. and J. F. Gerber. Estimating Evapotranspiration
from a Citrus Orchard with Weather Data, Proc Amer Soc for
Horticultural Science. Vol. 91, 1967.
171

69. Heaney, J. P. and W. C. Huber. Hydrologic Reconnaissance of
Conservation Areas I, 2 and 3 of the Central and Southern
Florida Project. Final report to Central and Southern Florida
Flood Control District, Dept. of Environmental Engineering,
University of Florida, Gainesville. 1971.
70. Heaney, J. P., A. I. Perez and J. L. Fox. Revised Nutrient Budget
Within the Organic Soils Area North of Lake Apopka. Report to,
East Central Florida Regional Planning Council, Dept. of
Environmental Engineering, University of Florida, Gainesville.
March 1972.
71. Henderson, F. M. Open Channel Flow. New York, The MacMillan Co.
1966.
72. Hills, R. C. "The Influence of Land Management and Soil Character
istics of Infiltration and the Occurrence of Overland Flow.
Journal of Hydrology. .13(2): 163181, June 1971.
73. Holtan, H. N. A Concept for Infiltration Estimates in Watershed
Engineering, U. S. Dept. of Agric., Agricultural Research Service
ARS 4151, 1961
74. Holtan, H. N. and D. E. Overton. StorageFlow Hysteresis in Hydro
graph Synthesis. Journal of Hydrology. £:309323, 1964.
75. Holtan, H. N. A Model for Computing Watershed Retention from Soil
Parameters. J of Soil and Water Conservation. 20£5):9l94, 1965.
76. Holtan, H. N., _et al_. Moisture Tension Data for Selected Soils
on Experimental Watersheds. ARS 41144, USDA Agric. Research
Service, Oct. 1968.
77. Hoopes, J. A. and D. R. F. Harleman. Waste Water Recharge and
Dispersion in Porous Media. MIT Hydrodynamics Lab Report No. 75
Cambridge, Massachusetts. June 1965.
78. Hoopes, J. A. and D. R. F. Harleman. Wastewater Recharge .and
Dispersion in Porous Media. Journal of the Hydraulics Division,
ASCE. 93CHY5):5171, Sept. 1967.
79. Hortenstine, C. C. and R. B. Forbes.. Soils Department, Institute
of Food and Agricultural Sciences (IFAS), University of Florida,
Gainesville, unpublished data, 1971.
172

80. Horton, R. E. The Interpretation of Runoff Plot Experiments with
Reference to Soil Erosion Problems. Proc Soil Sci Soc Am. 3^:340
349, 1938.
81. Horton, R. E. Approach Toward a Physical Interpretation of
Infiltration Capacity. Proc Soil Science Soc Am. £: 399417,
1940.
82. Hubbert, M.K. The Theory of Groundwater Motion. Journal of
Geology. 48:785944, 1940.
83. Huber, W. C. and A. I. Perez. Prediction of Solar and Atmospheric
Radiation for Energy Budget Studies of Lakes and Streams.
Publication No. 10, Water Resources Research Center, University
of Florida, Gainesville. 1970.
84. Huggins, L. F. and E. J. Moke. A Mathematical Model for Simu
lating the Hydrologic Response of a Watershed. Water Resources
Research. 4:529540, June 1968.
85. Hutchinson, G. E. Eutrophication, Past and Present. In: Eutro
phication; Causes. Consequences and Correctives. National
Acadamy of Sciences, Washington, D. C. 1969.
86. Jackson, D. R. and G. Aron. Parameter Estimation in Hydrology:
The State of the Art. Water Resources Bulletin. 7(3):457472,
June, 1971.
87. James, L. D., et^al.. An Evaluation of Relationships Between
Streamflow Patterns and Watershed Characteristics Through the
Use of OPSET  A Self Calibrating Version of the Stanford Water
shed Model. Research Report No. 36, University of Kentucky,
Water Resources Institute, Lexington, Kentucky. 1970.
88. Kaufman, W. J., et al. Disposal of Radioactive Wastes into
Deep Geologic Formations. Journal of the Water Pollution Control
Federation. 33(1):7384, Jan. 1961.
89, King, H.W. and E. F. Brater. Handbook of Hydraulics. New York,
McGrawHill, 1963.
90. Kisiel, C. C. Efficiency of Parameter and State Estimation
Methods in Relation to Models of Lumped and Distributed Hydro
logic Systems. (Presented at U. S.Japan Bilateral Seminar in
Hydrology. Honolulu, Hawaii. Jan. 1971.)
173

91. Kleinecke, D. Use of Linear Programming for Estimating Geo
hydrologic Parameters of Groundwater Basins. Water Resources
Research. 7(2) :367374, April 1971.
92. Klute, A. A Numerical Method for Solving the Flow Equation
for Water in Unsaturated Materials. Soil Science. 73:105
116, 1952.
93. Kohler, M.A. Lake and Pan Evaporation. In: Water Loss In
vestigations  Lake Hefner Studies. U.S.G.S. Prof. Paper No.
269. 1954.
94. Kunckle, S. Concentrations and Cycles of Bacterial Indicators
in Farm Surface Runoff. In: Relationship of Agriculture to
Soil and Water Pollution. Proc. Cornell University Cong, on
Agricultural Waste Management, Rochester, New York. Jan 1970.
p. 4960.
95. Kuznetsov, S. I. Recent Studies on the Role of Microorganisms
in the Cycling of Substances in Lakes. Limnology and Oceano
graphy. L3 (2):211224, 1968.
96. Li, W. H. and S. H. Lam. Principles of Fluid Mechanics. Reading,
Massachusetts, AddisonWesley Publ'. Co. May 1964. p. 5758.
97. Lichtler, W. F. and B. F. Joyner. Availability of Ground Water
in Orange County, Florida. Map Series No. 21. Florida Board
of Conservation, Tallahassee, Florida. 1966.
98. Lichtler, W. F., W. Anderson, and B. F. Joyner. Water Resources
of Orange County, Florida. Report of Investigation No^ 50,
U.S. Geological Survey, Tallahassee, Florida. 1966.
99. Linsley, R. K. and J. B. Franzini. WaterResources Engineering,
New York, McGrawHill. 1964.
100. Mackereth, F. J. Chemical Investigation of Lake Sediments
and their Interpretation. Proc Royal Society (B). 161:295
309, 1965.
101. Maddaus, W. 0. and P. S. Eagleson. A Distributed Linear
Representation of Surface Runoff. MIT Hydrodynamics Lab Report
No. 115, Cambridge, Massachusetts. June, 1969.
174

102. Matalas, N. C. Optimum Gaging Station Location. Proc IBM
Scientific Computing Symposium— Water and Air Resource Manage
ment, Yorktown Heights, New York. 1967. p. 8594.
103. Maxey, G.B. and J. E. Hackett. Application of Geohydrologic
Concepts in Geology. Journal of Hydrology. 1:3545, 1963.
104. McGauhey, P. H. and J. H. Winneberger. Final Report on a Study
of Methods of Preventing a Failure of SepticTank Percolation
Systems, SERL Report No. 6517, San. Engineering Research
Lab, University, of California, Berkeley. Oct. 31, 1965.
105. McGauhey, P. H., R. B. Krone and J. H. Winneberger. Soil Mantle
as a Wastewater Treatment System. SERL Report No. 667, San.
Engneering Research Lab, University of California, Berkeley,
Sept. 1966.
106. Miller, R. T., J. W. Biggar and D. R. Nielsen. Chloride Displace
ment in Panoche Clay Loam in Relation to Water Movement and
Distribution. Water Resources Research. J^: 6373, 1965.
107. Moore, S. F. The Design of Data Programs for Aquatic Ecosys
tems. Meeting Preprint 1574, ASCE Natl. Water Res. Engr.
Conf., Atlanta, Georgia. Jan. 1972.
108. Nakano, Y. and R. P. Murrman. A Statistical Method for Analysis
of Diffusion in Soils. Soil Sci of Amer Proc. 35(3):397402,
1971.
109. Nash, J.E. The Form of the Instantaneous Unit Hydrograph.
Inter Assoc Sci Hydro 1 Publ 45. 3_: 114121, 1957.
110. Nelson, R. W. InPlace Measurement of Permeability in Hetero
geneous Media. I. Theory of a Proposed Method. Journal of
Geophysical Research. 65j 17531758, 1960.
Hi. Nelson, D. W. and M. J. M. Romkens. Transport of Phosphorous
in Surface Runoff. In: Relationship of Agriculture to Soil
and Water Pollution. Proc. Cornell University Conf. on Agri
cultural Waste Management, Rochester, New York. Jan. 1970.
p. 215225.
112. Norton, T. E. and R. W. Hansen. Cattle Feedlot Water Quality
Hydrology. In; Animal Waste Management. Proc. Cornell Univer
sity Conf. on Agricultural Waste Management, Syracuse, New York
Jan. 1969. p. 203216.
175

113. Neuman, S. P. and P. A. Witherspoon. Variational Principles
for confined and Unconfined Flow of Groundwater. Water Resources
Research. £(5):13761382, Oct. 1970.
114. Odum, E. P. The Strategy of Ecosystem Development. Science.
164_:262270, 1969.
115. Onstad, C. A. and D. L. Brakensiek. Watershed Simulation by
Stream Path Analogy. Water Resources Research. £(5):965972,
Oct. 1968.
116. Orsborn, J. F. The Prediction of Piezometric Groundwater Levels
in Observation Wells Based on Prior Occurrences. Water Resources
Research. 2^1):139140, 1966.
117. Oster, J. D. and B. L. McNeal. Computation of Soil Solution
Composition Variation with Water Content for Desaturated Soils.
Soil Sci Soc Am Proc. 35_(3):436442, 1971.
118. Overton, D. E. Mathematical Refinement of an Infiltration
Equation for Watershed Engineering, ARS 4199, USDA Agricul
tural Research Service. 1964.
119. Overton, D. E. A LeastSquares Hydrograph Analysis of Complex
Storms on Small Agricultural Watersheds. Water Resources
Research. 4_(5): 955964, Oct 1968.
120. Park, R. A. and J. W. Wilkerson. Lake George Modeling Philoso
phy. Memo Report No. 7119. Rensselaer Polytechnia Inst., Troy,
New York, 1971.
121. Parizek, R. R., et^ al^. Waste Water Renovation and Conservation.
The Pennsylvania State University Studies, No. 23, University
Park, Pennsylvania. 1967.
122. Partheniades, E. Erosion and Deposition of Cohesive Soils.
Journal of the Hydraulics Division, ASCE. JU(HYI):105138,
Jan. 1965.
123. Peaceman, D. W. and H. H. Rachford. The Numerical Solution of
Parabolic and Elliptic Differential Equations. Journal of the
Society of Industrial and Applied Mathematics. 3(1):2841, March
1955.
176

124. Penman, H. L., Natural Evaporation from Open Water, Bare Soil and
Grass. Proc Royal Soc London, Series A. 1931120146, 1948.
125. Perez, A. I. A Water Quality Model for a Conjunctive Surface
Groundwater System. Doctoral Dissertation, Department of Environ
mental Engineering Sciences, University of Florida, Gainesville.
1972.
126. Philip, J. R., An Infiltration Equation with Physical Signifi
cance. Soil Science. 77_(2) : 153157, 1955.
127. Phillip, J. R. The Theory of Infiltration. I. The Infiltra
tion Equation and its Significance. Soil Science. 83:345
357, 1957.
128. Philip, J. R. The Theory of Infiltration. 4. Sorptivity and
Algebraic Infiltration Equations. Soil Science. £4:257264, 1957.
129. Putnam, H.D., £t aJL Eutrophication Factors in North Central
Florida Lakes. Publ. No. 5, University of Florida Water Resources
Research Center, Gainesville. 1969.
130. Remson, I., A. A. Fungaroli and G. M. Hornberger. Numerical
Analysis of Soil Moisture Systems. J Irrig and Drainage Div,
ASCE. 93(IR3):153166, Sept. 1967.
131. Remson, I., G. M. Hornberger and F. J. Molz. Numerical Methods
in Subsurface Hydrology. Wiley Interscience, New York. 1971.
132. Scalf, M. R., £t al_. Movement of DDT and Nitrates During Ground
Water Recharge. Water Resources Research. _5_(5): 10411052, Oct.
1969.
133. Scheidegger, A. E. Statistical Hydrodynamics in Porous Media.
J. of Applied Physics. 25_(8): 9941001, Aug. 1954.
134. Scheidegger, A. E. General Theory of Dispersion in Porous
Media. Journal of Geophysical Research. 66J10):32733278,
Oct. 1961.
135. Schneider, R. F. and J. A. Little. Characterization of Bottom
Sediments and Selected Nitrogen and Phosphorous Sources in
Lake Apopka, Florida. EPA Southeast Water Laboratory, Athens,
Georgia. March, 1968.
177

136. Schreiber, D. L. and D. L. Bender. Obtaining Overland Flow
Resistance by Optimization. Journal of the Hydraulics Division,
ASCE 98(HY5);429446, March, 1972.
137. Schultz, D. A. A Balance Sheet Method of Determining the Con
tribution of Agricultural Wastes to Surface Water Pollution. In:
Relationship of Agriculture to Soil and Water Pollution. Proc.
Cornell University Conf. on Agricultural Waste Management, Roches
ter, New York. Jan. 1970. p. 251262.
138. Sechrist, A., et^ al_. Mathematical Management Model—Unconfined
Aquifer. Report to Office of Water Resources Research, Water
Resources Center, Texas Tech University, Lubbock, Texas. 1971.
139. Seginer, I. A Model for Surface Drainage of Cultivated Fields.
Journal of Hydrology. 13(2):139152, June 1971.
140. Shamir, U.Y. and D. R. F. Harleman. Numerical and Analytical
Solutions of Dispersion Problems in Homogeneous and Layered
Aquifers. MIT Hydrodynamics Lab Report No. 89, Cambridge,
Mass. 1966.
i
141. Sheffield, C. W. and W. H. Kuhrt. Lake Apopka Its Decline and
Proposed Restoration. In: Proc of the Florida Env Eng Conf
on Water Pollution Control, University of Florida, Gainesville.
1969.
142. Sheffield, C. W. Effects of Agricultural Drainage Wastes.
(Presented at the Irrigation and Drainage Specialty Conference.
Miami Beach, Florida. Nov. 5, 1970.)
143. Simpson, E. S., C. C. Kisiel and L. Duckstein. SpaceTime
Samplings of Pollutants in Aquifers. (Presented at Symposium
on GroundWater Pollution, 15th General Assembly of Int. Union of
Geodesy and Geophysics, Moscow, USSR. Aug. 10, 1971.)
144. Sinha, L. K. and L. F. Lindahl. An Operational Watershed Model:
General Considerations, Purposes and Progress. Proc of Annual
Meeting of Am Soc Agric Eng (1970), Minneapolis Minnesota.
July 1970.
145. Smith, R, S. and D. A. Woolhiser. Overland Flow on an Infiltra
ting Surface. Water Resources Research. 7^4"):899913, Aug. 1971.
146. Soil Conservation Service (U.S.D.A.). Soil Survey—Orange County,
Florida. Gainesville, Florida. Sept. 1960.
178

147. Soil Conservation Service (U.S.D.A.). General Soil Map—Lake
County, Florida and General Soil Map—Orange County, Florida.
Gainesville, Florida. May 1967.
148, Stanford, G., C. B, England and A. W. Taylor. Fertilizer Use
and Water Quality. USDA Agricultural Research Service, ARS
41168. Oct. 1970.
149. Stephens, J. C. Peat and Muck Drainage Problems. J Irrig and
Drainage Div, ASCE. S£(IR2):, June 1969.
150. Stemberg, Y. M. Parameter Estimation for Aquifer Evaluation.
Water Resources Bulletin. 7^3):447456, June 1971.
151, Stewart, E. H., D. P. Powell and L. C. Hammond. Moisture Charac
teristics of Some Representative Soils of Florida. ARS 4163.
USDA Agric. Research Service, April 1963.
152. Stewart, E. H., et_ aJL Evaporation as Affected by Plant Density
and WaterTable Depth. Am Soc Agric Engr Trans. 12(5):646647,
1969. 
153. . Stone, H. L. and P. L. T. Brian. Numerical Solution of Convection
Transport Problems. Journal American Inst of Chemical Engineers.
£(5):681688, Sept. 196T~
154. Stratton, F. E. Nitrogen Losses from Alkaline Water Impoundments.
Jour San Eng Div, ASCE. £5(SA2):223231, April 1969.
155. Tanji, K. K., e£ a^. II. A Computer Method for Predicting Salt
Concentrations in Soils at Variable Moisture Contents. Hilgardia.
J58(9):307318, June 1967.
156. Theis, C. V. The Relation Between the Lowering of the Piezo
metric Surface and the Rate and Duration of Discharge of a Well
Using GroundWater Storage. Trans Amer Geophysical Union. 16:
519524, 1935.
157. Thiem, G. Hydrologische Methoden. Leipzig, Gebhardt. 1906.
158. Thornthwaite, C. K. and J. R. Mather. The Water Balance.
Laboratory of Climatology Report, Drexel Inst. of Technology,
Centerton, New Jersey. 1955.
179

159. Thornthwaite, C. W. and J. R. Mather. Instructions and Tables
for Computing Potential Evapotranspiration and the Water Balance.
Lab. of Climatology Report, Drexel Inst. of Technology, Centerton,
New Jersey. 1957.
160. Thomas, L. H. Elliptic Problems in Linear Difference Equations
over a Network. Report of Watson Science Computing Lab, Columbia
University, New York. 1949.
161. Todd, D. K. Ground Water Hydrology, 6th ed. John Wiley and Sons.
Jan. 1967.
162. Toth, J. A Theory of Ground Water Motion in Samll Basins in
Central Alberta. Journal of Geophysical Research. 67(11):
43754387, 1962.
163. Toth, J. A Theoretical Analysis of Ground Water Flow in Small
Drainage Basins. Journal of Geophysical Research. 68_(16):
47954812, 1963.
164. Townshend, A. R., K. A. Reichert and J. H. Nodwell. Status
Report on Water Pollution Control Facilities for Farm Animal
Wastes in the Province of Ontario. In: Animal Waste Management.
Prbc. Cornell University Conf. on Agricultural Waste Management,
Syracuse, New York. Jan. 1969. p. 131149.
165. Trelease, F. J, and M. W. Bittinger, Mechanics of a Mathmatical
Groundwater Model. J Irrig and Drainage Div, ASCE. 89(l):5162,
March 1963.
166. U. S. Army Corps of Engineers. Everglades Gaging Program Pro
gress Report—Everglades Area, Florida. Report No. 6. U. S.
Army Engineer Dist., Meteorology and Regulation Section, Jackson
ville, Florida. Dec. 1956.
\
167. U. S. Dept. of Agric. and U. S. Dept. of Commerce. Monthly
Precipitation Probabilities by Climatic Divisions. Miscellaneous
Publ. No. 1160. 1969.
168. U. S. Weather Bureau. Rainfall Frequency Atlas of the U. S.
Technical Paper No. 40. May 1961.
169. Warrick, W. A. J. W. Biggar and D. R. Nielsen. Simultaneous
Solute and Water Transfer for an Unsaturated Soil. Water Resources
Research. 1(5):12161225, Oct. 1971.
180

170. Weaver, H. A. and W. H. Speir. Applying Basic Soil Water Data
to Water Control Problems in Everglades Peaty Muck. ARS 4140,
Agric. Res. Service, U.S.D.A. Nov. 1960.
171. Webber, L. R. and T. H. Lane. The Nitrogen Problem in the Land
Disposal of Liquid Manure. In: Animal Waste Management. Proc.
Cornell University Conf. on Agricultural Waste Management, Syra
cuse, New York. Jan. 1969. p. 124130.
172. Weibel. S. R., e£ al.. Pesticides and Other Contaminants in Rain
fall and Runoff. J. AWWA. 518:10751084, 1966.
173. Williford, J. W., e£ al_. The Movement of Nitrogenous Fertilizers
through Soil Columns. In: Collected Papers Regarding Nitrates
in Agricultural Waste Waters. FWQA, Dec. 1969.
174, Winterer, E. V. Percolation of Water Through Soils. Calif.
Aqric. Exp. Sta. Annual Rept. from July 1, 1920 to June 20,
1922.
175. Yeh, W. and G. W. Tauxe. Quasilinearization and the Identification
of Aquifer Parameters. Water Resources Research. 2.C2) :375381,
April 1971.
181

SECTION XI
PUBLICATIONS
Perez, A. I., W. C. Huber, J. P. Heaney and E. E. Pyatt. A Water
Quality Model for a Conjunctive Surface— Groundwater System.
(Presented at the Seventh American Water Resources Conference,
Washington, D. C. October 1971.)
Perez, A. I., W. C. Huber, J. P. Heaney and E. E. Pyatt. Back
ground of a Surface— Groundwater Quality Model. (Presented at ASCE
Water Resources Engineering Conference, Atlanta, Georgia. January
1972.)
Perez, A. I. A Water Quality Model for a Conjunctive Surface—
Groundwater System. Doctoral Dissertation, Environmental Engineering
Sciences Department, University of Florida, Gainesville. June 1972.
Perez, A. I., W. C. Huber, J. P. Heaney and E. E. Pyatt. A Water
Quality Model for a Conjunctive Surface— Groundwater System: An
Overview. Water Resources Bulletin. 8^(5): 900908, October 1972.
182

SECTION XII
LIST OF SYMBOLS
a = constant in infiltration equation; also, constant in
dispersion equation
A = infiltration parameter
A, B, C = lefthand side coefficients in tridiagonal matrix
AC, BC, CC = lefthand side coefficients in tridiagonal matrix
c • concentration; also, a subscript indicating corrected
value
c* = intermediate value of concentration
C B slope of moisture versus pressure head diagram
C B intrinsic permeability
d B soil particle diameter
d = thickness of "equivalent layer" for mole drains
o
e
D B righthand side coefficient in tridiagonal matrix; also,
dispersion coefficient
DC = righthand side coefficient in tridiagonal matrix
D = distance
5
D . D = dispersion coefficients in x and ydirections
x y
e = vapor pressure
E = evaporation rate
f = infiltration rate
f. = initial infiltration rate
f = final infiltration rate
F = potential volume of infiltration
P
g = gravity
183

h = distance between water table and drain
i = rainfall intensity; also, cumulative infiltration
i, I = column index
i = rainfall excess, (rainfall less infiltration)
I = total infiltration; also, infiltration rate
j, J  row index
k = Shield's parameter; also, decay coefficient
K = hydraulic conductivity or permeability
K_ = dispersion coefficient
L ~ catchment length; also, length of drains
m = exponent constant in kinematic wave formulation
(=5/3 for turbulent flow)
n = time step; also, Manning's roughness coefficient; also,
a constant; also, a subscript indicating new value
p = subscript indicating predicted value
P = mass of pollutant; also phosphorus
P = initial mass of pollutant
q = overland flow per unit width
q = flux of sediment
S
Q = flow rate; also groundwater recharge
r  nutrient release rate
R = rainfall rate; also rainfall less evaporation; also,
hydraulic radius
R = total rainfall
R = distance ratio
s
S = source/sink; also, sorptivity; also mole drain spacing
S, = bed slope
D
S. = sink rate of nutrients from water
i
S = specific gravity
S_ = source rate of nutrients to water
t = time
t = time of concentration
c
184

tf = duration of infiltration
t = rainfall duration
u = nutrient uptake rate
u, v, w = velocities
v = critical scouring velocity
c
w = relaxation factor in groundwater flow calculation; also,
water withdrawal rate per unit depth
W = water withdrawal rate; also, weight
x, y, z  space coordinates
y = flow depth in overland flow
y = flow depth at catchment outlet
Li
z = elevation head
a = square of ratio of horizontal to vertical spacing in
groundwater model; also, constant in kinematic wave over
land flow computations; also, decay constant in infiltra
tion equation
a_, 3 = constants in dispersion equation
•5 3
Y = unit weight of water
6 = volumetric moisture content (fraction); also angle of
catchment slope
v = kinematic viscosity of water
p = density of water
= hydraulic head
ty = pressure head; also, constant in sediment load calculations
T = shear stress at the stream bed
T = critical shear stress for erosion
cr
A list of computer variable names is provided in the comments section
at the beginning of each program.
185

SECTION XIII
APPENDICES
Page
A. Details of Surface Water Models 187
B, Details of Unsaturated Zone Model 201
C. Details of Groundwater Model 225
D. Details of Interfacing Methodology 232
E. Details on Application to the Lake Apopka Basin 236
F. Summary of Additional References on Agricultural Nutrients 257
G. Listing of Computer Programs 259
186

APPENDIX A
DETAILS OF SURFACE WATER MODELS
INTRODUCTION
Since adequate information is available on the Stormwater Management
Model runoff and receiving portions (Environmental Protection Agency,
1971), the discussion on surface water models is confined to the kine
matic wave overland flow formulation as described in Section IV. The
mathematics of this formulation are covered thoroughly by Eagleson
(1970) and, therefore, the emphasis herein is on illustrations and
sample calculations. Quality computations are included. Finally, a
description of the various program subroutines is provided. Consistent
English units are used throughout. Metric equivalents are not given
since the data are hypothetical.
SAMPLE CALCULATIONS
Catchment Geometry
The geometric characteristics of a hypothetical catchment are illus
trated in Figure 73. Part (a) of this figure is a plan view of the
catchment and its boundaries and assumed dimensions are indicated.
A side view of the catchment is presented in part (b) on which the
slope is shown.
Hydraulic Parameters
As discussed in Section IV, the flow rate at the catchment outlet can
be computed from the formula:
q  °tym (42)
where q = discharge per unit width, ft3/sec/ft
y = depth of water, ft
a and m = constants
187

oo
oo
TOPOGRAPHIC
CONTOUR
LINES
SURRVCE BODY
OF WATER
PART A  PLAN VIEW
f f FTRAINFALL
100 ft
PART 8  LATERAL VIEW
Figure 73. Plan and lateral views of a hypothetical catchment.

The values of a and m vary according to whether the flow is laminar or
turbulent. Eagleson has pointed out that since both flow regimes are
likely to occur on vegetated surfaces, a compromise value of m « 2
might be used (Morton, 1938) . The corresponding value of a must then
be determined from field measurements,
The lack of appropriate data in the Lake Apopka basin, and possibly many
other areas, precludes the statistical estimation of the parameter a.
For this reason, the flow was assumed to be turbulent (m = 5/3) . This
assumption permitted the direct calculation of from:
(43)
where 6 = catchment slope angle
n = Manning roughness coefficient
In this case, the units of a are (ft1/ '/sec).
The value of n can be estimated from tabulated values for different
types of land surfaces (Henderson, 1966). In vegetated areas, it appears
that n varies with flow depth (U. S. Army Corps of Engineers, 1956). For
simplicity, a fixed value of n • 0,05 is assumed for the hypothetical
catchment in question.
Hydrologic Input Conditions
The hydrologic inputs required to perform the overland flow calculations
are rainfall rates, infiltration rates and ponding depths. The latter
two are obtained from the unsaturated zone model. The assumed values
of these three variables are indicated in Figure 74. Evaporation rates
are not considered.
Preliminary Calculations
Since the analytical kinematic wave formulation requires the usage of
constant rainfall and infiltration inputs, preliminary calculations
must be conducted to obtain timeaverage values for these variables.
The first step is to inspect the ponding depth values to determine the
time at which maximum depth occurred. This time is then assumed to be
the duration, tr, of the simplified rainstorm, which is 15 min (Figure
74).
Next, the total rainfall, R, is computed as follows:
189

LJ
CC
Sc
e
7
i
6
5
3
2
£
\


•
1
1
t 1 1 III
10
IS
20
25
30
< c
(E ""_
£ CK
O
4
3
1
(

] 1
i i i i ili
35 10 15 20 25 30
0. ?
UJ 2
1  1
10
1 l
15
20
25
30
TIME, min
Figure 74. Hypothetical rainfall rates, infiltration rates
and ponding depths used in overland flow
calculations.
190

inmin
[8(2.5} + 7(5) * 4(5) + 2(5) + 1(5) * 0.5(5)] hr _ .
60 min/hr  lt54 in
and the synthetic rainfall intensity becomes
= 0.102 in/rain = 0.00858 ft/min f or 0 < t < t
— — r
= 0 for t > t
Turning to the infiltration calculations, it is assumed that the time
at which ponding disappears constitutes the end of the infiltration
process. An inspection of Figure 74 reveals that the duration of this
process is tf » 27.5 min.
The total amount of infiltration for that period is
inmin
j m [1(2.5) + 1.5(5) * 1.5(5) + 2(5} + 0.5(5) + 0.5(5)] hr
T "" 60 min/hr
* 0.538 in
Accordingly, the synthetic infiltration rates becomes
= 0.0196 in/min = 0.00162 ft/min for 0 £ t 5 tf
for t > t£
Kinematic Wave Calculations
Two main cases arise in the formulation, depending on the relative
values of the storm duration, tr, and the time of concentration, tc.
The latter is computed from the formula:
191

Li
1/m
(44)
where L = catchment length, ft
i^ = rainfall excess (rainfall less infiltration in ft/min)
Substitution of the appropriate parameters in (44) yields:
c
100(0.00858
1.49flO V
0.05[100 J
 0.00162)'2/3
1/2 f60 sec)
[ minj
3/3 f ,
(2800)
[565J
0.6
=2.6 min
The correction factor in the denominator was necessary to convert a from
(ft1/3/sec) to (ft1'3/min). For the particular catchment and hydrologic
conditions in question, t exceeds t .
JL C
It is now possible to compute the depths of water, and the correspond
ing velocities and flows, at various times. These values are computed
herein only at the catchment outlet. The calculations for the first
time step (t = 5 min) are given below.
The equations used to obtain the flow depth are
y = (i  f)t 0 £ t £ t
(45)
(46)
A more complex relationship must be used for the case t >tr. Since in
this case t < t < t ,
£ — — f
the depth becomes
(0.00696)(2.6) = 0.0181 ft
192

The corresponding velocity is
v = ay51"1 = (9.42)(0.0181)°>66 = 0.705 ft/sec
and the flow per unit width becomes
q = aym = vy = 0.01275 ft3/sec/ft.
Multiplying this unit flow by the width (ft), the resulting catchment
outflow is
Q = qw = (0.01275)(60) = 0.765 ft3/sec.
Sediment Load Calculations
The sediment transport formula used in this study is that developed by
duBoys (Einstein, 1964):
To Tcr
Y (47>
' t
where q = flux of sediment, Ib/sec/ft of width
ty = a constant for a given bed composition
T a bed shear stress, lb/ft2
o
T a critical shear stress, lb/ft
cr
2
3
Y a unitweight of water, lb/ftj
The bed shear stress can, in turn, be expressed as:
T = yys (48)
where y = flow depth, ft
S, = slope of the bed
193

Substitution of equation (48) into (47) yields
q
s
T
Y
(49)
which is the actual expression utilized herein. Values of if> and Tcr
are tabulated (Einstein, 1964) for various grain sizes.
Assuming that the soil in the hypothetical catchment is fine sand,
the sediment load at the first time step (t = 5 min) is:
n — fc o*z v i n5> /*i QI v
n * ID.^J ^ J.U ) I J. • o J. ^
Te V** • **** *fcv ^ V ^ • w * — — ,/ V ~ "" ^1***** •* v x « 4
s I 62.4
=1.47 Ib/sec/ft
Thus, the total catchment load becomes
Qs = wqg =88.3 Ib/sec.
The calculations on the nitrogen and phosphorous fractions in the
sediment load are shown in Appendix E.
PROGRAM DESCRIPTION
The kinematic wave overland flow program is composed of a main program
(MAIN) and several subroutines, AVE, WAVE, RECV, HIDUR, LODUR, RQUAL
and WRIT. All these subroutines are involved in quantity computations,
with the exception of RQUAL and WRIT. The former handles water quality
computations and the latter writes out the overall results.
Main Program
The first task of MAIN is to read in the necessary hydrologic, geometric
and quality data. Included among the geometric information is the
numbering system that connects the catchments with the receiving elements.
Following the data input, MAIN calls subroutines AVE, WAVE^ RECV, RQUAL
and WRIT, in that order. Subroutine AVE performs the averaging compu
tations for rainfall, evaporation and infiltration, and thus prepares
these data for use in the kinematic wave calculations. The task of
WAVE is to carry out these calculations for all catchments and time
steps. It ultimately computes flows into the receiving elements.
Subroutine RECV then computes the total flows into each receiving element
194

from all its feeder catchments. The functions of RQUAL and WRIT were
explained above.
Subroutine AVE
This subroutine determines the mean values of infiltration and rain
fall for the time period considered. The stepwise procedure used to
obtain these values, discussed below, are also shown in Figure 75 and 76.
!
Beginning with infiltration, AVE first determines the time step after
which this process no longer occurs. The method used is to scan the
infiltration results from the unsaturated zone program and establish
when the ponding depth becomes zero. The time interval corresponding
to this event is considered to constitute the end of the infiltration
process. The scanning process described above is illustrated as Step 1
in Figure 75.
In actual practice the runoff process will probably produce a zero
ponding depth condition at an earlier time than that obtained from the
above scanning method. The reason is that the unsaturated zone model
assumes a laterally stagnant, nondynamic, ponding situation. Never
theless, it is thought that the approximation used is adequate for the
modeling of rural watersheds. In the future, one method to obtain
more accurate results might be to interface the infiltration and over
land flow models at every time step.
Following determination of the duration of infiltration, AVE computes
the timeaverage value of the infiltration rate. This calculation is
depicted as Step 2 in Figure 75. Finally, updated values of the infil
tration rate are stored for all the time intervals in the study period,
in the manner shown in Step 3 of the same figure.
Subroutine AVE then considers the rainfall process, and converts the
original irregularshaped hyetograph into a step function containing
the same total amount of rainfall. This conversion process is begun
by scanning the output from the unsaturated zone program, determining
the time at which the maximum ponding depth occurs and arbitrarily
setting the duration of the synthetic storm equal to the above time.
This process is illustrated in Step 1 of Figure 76. The subroutine
then computes the total rainfall from the original record, as shown
in Step 2 of the same figure. Finally, it constructs a synthetic step
function hyetograph (Step 3) using the duration, tr, determined above.
Subroutine WAVE
To better represent the overland flow process, it is desirable to apply
the .kinematic wave relationships not only at each catchment outlet,
but also at several points along the length of the catchment. For this
195

STEP I SCAN PONDING DEPTHlWATER DEPTH)
DATA FROM INFILT TO DETERMINE
END OF INFILTRATION PERIOD.
2
ASSUMED END OF
INFILTRATION PROCESS
23 V Ml M
TIME INTERVAL NUMBER
i
UJ
I
O
I
b
i
STEP 2 USE THE INFILTRATION RATES FOR M
PERIODS ABOVE FROM INFILT TO
OBTAIN AN AVERAGE VALUE.
2 3 V Ml
TIME INTERVAL NUMBER
M
STEP 3 USING THE RESULTS OF (2) ABOVE
CONSTRUCT AN EQUIVALENT STEP
FUNCTION.
2 3 ^l M
TIME INTERVAL NUMBER
M + l M+2
Figure 75. Graphical representation of computational steps
in subroutine AVE involving infiltration.
196

MAXIMUM DEPTH?
STEP I SCAN PONDING DEPTHS FROM INFILT
TO DETERMINE TIME OF MAXIMUM
PONDING DEPTH
f
DEPTH.
1 1
23 4 5 6 78
TIME INTERVAL NUMBER
I
SHADED AREA IS TOTAL RAINFALL
TEP2 USE ORIGINAL HYETOORAPH TO
COMPUTE TOTAL RAINFALL.
2 3 4 f 5 6 789
TIME INTERVAL NUMBER
UJ
*
a:
ce
SHADED AREA SAME AS
IN (2) ABOVE
STEP 3 CONSTRUCT SYNTHETIC RAINFALL IN
STEP FUNCTION FORM.
123456 7 •
TIME INTERVAL NUMBER
Figure 76. Graphical representation of computations
in subroutine AVE involving rainfall.
197

purpose each catchment is divided into a given number of intervals,
and the formulation is applied to the downstream end of each interval.
This division process is illustrated in Figure 77. In this study, only
one interval, comprising the entire catchment, was considered. This
was acceptable because of the small size of each catchment.
The type of formulas used at each catchment location largely depends
on the relative values of tr, the storm duration as defined in AVE,
and tc, the time of concentration. The latter term can be defined as
the maximum amount of time during which an increase of depth can occur
on the surface. While tr is assumed to be the same throughout the
basin, tc varies with location. At each time step and each catchment
interval, WAVE compares these values. If the storm has a relatively
high duration, namely tr exceeds tc, subroutine HIDUR is called. Con
versely, if tc is greater than tr, LODUR is utilized. Each of these
two subroutines computes the depth of water at each catchment location.
After the flow depths have been determined, WAVE proceeds to calculate
their corresponding velocities and flows by means of equations (7)
and (6). These values are then stored for future use.
Subroutine RECV
At each time step, this subroutine computes the total flow entering
each receiving element from its feeder catchments. To accomplish this
purpose, it considers the stored flows corresponding only to the down
stream interval in each catchment.
Subroutine HIDUR
The kinematic wave computations carried out in this subroutine corres
pond to the case tr > tc as mentioned above. These computations branch
out according to three possible cases. To determine which particular
case holds, HIDUR compares the relative values of t, the time corres
ponding to the time step in question, to three other parameters: tc, tr
and t^. The latter parameter is defined as the time at which the depth,
y, at a given location vanishes. At the conclusion of these calculations,
control is returned to subroutine WAVE.
Subroutine LODUR
This subroutine is utilized in a relatively low duration situation;
namely, tr < tc. The computations herein are broken down into three
classes, depending on the relative values of t^, defined above, and
tp, which represent the time at which the depth begins to decrease.
Each one of these three cases is in turn divided into three subcases.
As in HIDUR, control is returned to WAVE dollowing completion of the
computations.
198

RAINFALL
to
to
INFILTRATION,
RECEIVING
ELEMENT
Figure 77. Lateral view of overland flow situation, indicating catchment divisions.

Subroutine RQUAL
This subroutine utilizes the duBoys sediment transport equation (Einstein,
1964) to compute the suspended solids fluxes at the catchment outlets.
It also calculates the nitrogen and phosphorous fluxes at the outlets.
Subroutine WRIT
After subroutine RQUAL has completed its calculations, subroutine WRIT
is called to print out the results obtained. The most important re
sults are the flows and concentrations at the various catchments, and
the total inputs to the receiving elements.
2<)0

APPENDIX B
DETAILS OF UNSATURATED ZONE MODEL
INTRODUCTION
The first item to be discussed is the derivation of coefficients used
in the matrix solutions of both the flow and quality portions of the
unsaturated zone model as developed in Section V. Subsequently, sample
calculations are provided for one time step. Finally, a discussion of
the various subroutines is given.
DERIVATION OF COEFFICIENTS
Flow Model
The basic linknode structure used in a given soil profile is shown
in Figure 78. Also indicated therein are the boundary conditions at
the top and bottom of the profile.
The flow governing equation used in this study was that developed by
Freeze (1969), with an additional sink term:
s = J
3t 3z
(50)
where ip = pressure head, cm
C(ip) = slope of the moisture versus pressure head curve, I/cm
t = time, min
S = sink, 1/min
z = elevation, cm
K(i0 = hydraulic conductivity, cm/min
201

TOP BOUNDARY CONDITIONS:
K>
o
to
RAINFALL
EVAPORATION
GROUND
SURFACE
S////J///S//SSSSSSS S TX / X / /^//SSSSSS ' S S
Ll
I
NODE
NUMBERS.
.WATER TABLE
BOTTOM BOUNDARY CONDITION: RECHARGE
Figure 78. Illustration of unsaturated zone profile, showing nodes
and boundary conditions.

The sink terra above includes withdrawals by plants and artificial
drains. Its usage is illustrated in Appendix E.
jnterior Node—
The finitedifference expression corresponding to equation (50), written
for an interior node (1 < j < L), is:
;M
At
_L
Az
Az
n
n
(51)
in which the subscripts denote the node and the superscripts the time
step. The expressions for ^,, if)__, and tyrrr are the following:
+ *?*. *
(52)
*
II
,n
(53)
(54)
These expressions are valid only for n > 2. Otherwise, the definitions
Provided by Rubin and Steinhart (1963) should be used.
203

Since the flow governing equation form desired is:
V?1
it becomes necessary to regroup terms in equation (51). Expanding
first, the following equation is obtained:
C55)
2Az
2Az
2Az
2Az
2Az
2Az
2Az
^ (56)
Rearranging terms and multiplying both sides by (1), the expression
becomes
2Az
.. .
2Az 2Az
At
2Az
2Az 2Az
(57)
Note that the terms on the righthand side are known quantities.
Comparing equation (57) to (55) it follows that:
2Az
(58)
204

Bi = 2Az * 2Az* * At
n
2Az
(59)
(60)
n
At
2Az 2Az
(61)
The expressions for the four coefficients above are valid for any in
terior node. However, the expressions for the top and bottom nodes
constitute special cases.
Top Node
The boundary condition at the ground surface (j « L) is (Freeze, 1969):
JfWU
 1
(62)
where R « net rainfall infiltration (rainfall less evaporation),
cm/rain
If a sink term is included, equation (62) can be rewritten as:
R  S(Az/2)  K
(63)
Where • ijj + z
The finitedifference expression for (63), after some rearranging,
becomes
205

s"(Az/2)Az
K
r i )
n
' '"I ^ *L1 "
K
r i \
nT
Thus, the coefficient expressions become
A?  0
L ~ L
(64)
(65)
(66)
RAz
" ZL + ZL1
S"(Az/2)Az
Li
(67)
Bottom Node—
The boundary condition at the bottom node (j = 1) can be written as:
K
3z
= Q + s"(Az/2)
(68)
Where Q = groundwater recharge (positive downward), cm/min
After some manipulations, the following expression is obtained:
QAz . „ S1(Az/2)Az
K
, j »
n y
% .
• *1 '2
K
r i ^
nf
* i
(69)
206

It follows that the coefficient expressions are:
cn = o
QAz
""*
Zl + Z2
(70)
(71)
(72)
Quality Model
The quality governing equation consists of a convectivedispersion
formulation with provision for source/sink reactions:
where
KnCv)
S  S.
o i
D
(73)
time> min
moisture fraction, cm/cm
concentration, mg/1
velocity (positive upward), cm/min
elevation, cm
dispersion coefficient, cm2/min
source, mg/l/min
sink, mg/l/min
Since the flow through the soil will be downward following rainstorms,
it is convenient to redefine the velocity sign convention as positive
downward. The resulting expression becomes:
207

v * VCV) t S  S. (74)
9t 9z D ' 3z2 01
The computation of source/sink terms is illustrated in Appendix E. To
solve (74) numerically, it is necessary to regroup its terms and obtain
an expression of the form:
V?  (AB)ncn+i = (DC)J (75)
The expressions for coefficients AC, BC, CC and DC are derived below.
Interior Node—
For an interior node (1 < j < L), equation 74 can be written in finite
difference form as:
(en  eH n (cn  c
"1 3 3 J + en I J J
j At j At
n n , n n
'
After some algebra, the following expressions are obtained
vn
en _ 9ni en 2KD
BCn 
BC
.
j « it C78)
v"
cc" '  + C793
208

rn _ _ n
UL. = o
(80)
Top Node
At the uppermost node (j = L), the dispersion and source/sink terms are
neglected. Its concentration at the first time step is set, and its
value at subsequent steps is computed as follows:
1. Determine the distance, Ds, that the substance traveled between the
current and the previous time step.
(81)
2. Calculate the ratio, Rs, of this distance to the midpoint distance
between the two upper nodes
R =
s (Az/2)
(82)
3. Develop a massbalance relationship involving the current and pre
vious time steps:
L L
L L S
(83)
4. Substitute (82) into (83) and rearrange terms, to yield:
n
L ni
en L
vl^At
i ^
f Az/2 J
(84)
Thus, coefficients AC, BC, CC and DC are not used explicitly to compute
the concentrations at the top node.
209

Bottom node—
Due to the absence of two adjacent nodes, required to evaluate a second
derivative, the dispersion term is not considered at the bottom node
(j « 1). The finitedifference form of equation (74) at this location
is:
ffln ftni) f n nil f n n
niei  91 J + en K " c! J _ n [C2 " cl
1 At °1 At " vl I
C  Si
°i h
(85)
The resulting expressions for the coefficients, after regrouping of
terms,
(86)
At
At Az
(87)
C? = 0
(88)
c ni
 s
n
(89)
Discussion of Dispersion Equation—
The expression used to compute the dispersion coefficient was derived
by Hoopes and Harleman (1965) for sands based on laboratory experiments.
It is unlikely to be generally applicable, but was useful for this
study. The expression is:
(90)
V
where v = kinematic viscosity of water, craVsec
a_ = 88 (for sand)
«J
v = velocity, cm/sec
210

C = coefficient of intrinsic permeability, cm2
B = 1.2
The usage of this formula requires that the computed velocities at the
nodes, which have units of Ccm/min) be converted to (ft/sec).
Intrinsic permeabilities are computed from values of the hydraulic
conductivity. Since the kinematic viscosity of water is a function
of temperature, the value corresponding to the mean annual temperature
in the study area was used.
After KD is computed, it must be converted to units of cm2/min for
usage in the quality calculations. The various conversion factors used
in connection with equation (90) are given in the listing of the pro
gram (see Appendix G) .
SAMPLE CALCULATIONS
In both the flow and quality models, the calculations follow a twostep
procedure: 1. compute the values of the coefficients of each node, and
2. solve the resulting system of equations using Gaussian elimination
(Hadley, 1964). Step (1) is illustrated below for a single node in
the soil profile. Space limitations preclude the illustration of step
(2) . The soil type considered herein was the sand of the Lake Apopka
basin, whose estimated hydraulic characteristics are shown in Figures
42 and 86.
Flow Model
An isolated portion of a soil profile is shown in Figure 79. The
calculations to be discussed are those for node j = 5 and time step
n « 7. Indicated in the figure are the pressure heads for this node,
and its two adjoining nodes, at the previous time step (n = 6) . Since
the computations of withdrawals (sinks) by plants and artificial drains
are explained in Appendix E, these sinks are considered to be zero in
this discussion.
The coefficients AL,Bg, C75 and D* can be obtained from the formulas
for an interior node, that is, equations (77), (78), (79) and (80).
A review of these expressions indicates that to determine certain values
of K and C, it is first necessary to compute ty*, ^TT7 and ^TTT?
5 5 5
Equations (52), (53) and (54) for these three variables contain terms
for pressure heads at the present time step, which are yet unknown.
Thus these values must be estimated by any one of several methods,
211

to
»«
tsJ
TIME STEP
B=6
TIME STEP
„= 7
J JJ S s S J J S J
PRESSURE HEADS cm
ATn=6
V4 = 60.1
Y5 = 72.4
V6 = 80.3
COORDINATE AXES:
L
>J/J*V*>
.
'
GROUND SURFACE
•rs s s f f f f
j=5
J=4
o—
AzIOcm
1=
?
At * lOmln.
WATER TABLE
Figure 79. Illustration of an isolated portion of an unsaturated zone system, showing
flow model data.

including regression. For simplicity, these pressure heads are assumed
to be equal to those in the previous time step. Substitution into
equations (52), (53) and (54) yields:
cm
, 7 f 72.4  80.3  80.3  72.4) _. ,
^ =    . _ 76i3
5 '
, . [ ^O.l  72.44 72.4  60.1J . . w>2
5 *" '
7 (• 72.4  72.4)
"'ill B [  2  J a  72'4
The conductivities and slope parameters corresponding to .these inter
mediate head values can be estimated from Figures 86 and 42. The
resulting coefficients are:.
D;  (0.0038,( 80.3) * [(0.00378) (10). 0.076. 0.096)^ „.,,
+ CO. 0048) ( 60.1) + 0.076  0.096 + 0 •  0.266
Thus, the basic finitedifference equation for j • 5 and n • 7 is:
 (0.0048)l(i7 + (0.0124W7  (0.0038)i»7 = ( 0.266) (91)
it 5 6
Following the computation of pressure heads throughout the soil profile
using the elimination technique, it is possible to determine the flow
velocity at each node. The finite difference formula for Darcy's law,
applied to node j = 5, is:
213

C92)
where v denotes the velocity (cm/min) , which is defined as positive in
the downward direction, and z the elevation (cm) . It may be noted that
Z  Z = 2Az =  20.0 cm
For demonstration purposes, it may be assumed that $1 = 56.1, if;7.
 69.0 and i^7 =  78.1. From Figure 86, the conductivity becomes
K7 = 0.074 cm/min
Substitution into equation (92) produces
v' .  (0.074) [( 56a>  I" "•" * < = 0.0074 cm/min
5 ^
Quality Model
The computations of nutrient withdrawals by pjlants and drains, and of
source/sink reaction rates are shown in Appendix E. For this reason,
the terms involving these processes are set equal to zero below.
An isolated portion of a soil profile is depicted in Figure 80. Shown
therein are assumed values of moisture content, 9, for three nodes at
two time steps. Also indicated is the velocity, v, at the middle node
at the current time step n = 10, and the concentration at n = 9. In
practice, these moisture and velocity values are obtained from the flow
model. The concentrations at the various nodes at the current time
step constitute the unknowns.
To compute the values of the coefficients AC, BC, CC and DC at node
j =4, it is first necessary to determine the dispersion coefficients
at that location. As can be seen in equation (90), the computation of
this coefficient requires the values of the kinematic viscosity and
hydraulic conductivity.
Values for the kinematic viscosity of water are available for different
temperatures (King and Brater, 1963). The subsurface water temperature
214

ro
t>
tn
TIME STEP TIME STEP
n = 9 n=IO
GROUND
/ / / s / /
VALUES
®!
e
el
c!
VALUES
v,"
®4
.2
_
// /"/ / S / / S S S / /
 0.24
= 0.20
= 0.16
= 5.2 ma/I i  R
ATn=lO= J = 4
= 0.06 cm/mi n . _ 3
= 0.001 cm/iec
= 0.25
mnornuATc
CUUKUI NA 1 1
AXIS
t
1
<:
b
<
, i
^ ^
L f
<
At = lOmin.
>
>
r
Az = 10 cm
?
/
/WATER
TABLE
Figure 80. Illustration of an isolated portion of an unsaturated zone system,
showing quality model data.

in the Lake Apopka basin was estimated to be 71.5°F, which is the
annual mean value for Orlando, Florida (Lichtler, Anderson and Joyner,
1968). Accordingly,
v = 1.040 x io~5 ft2/sec = 1.12 x 10~8 cm2/sec
To determine the hydraulic conductivity, it is necessary first to
determine the pressure head corresponding to the moisture content at
j = 4. From Figure 42,
=  59.0 cm
The corresponding hydraulic conductivity value, K, obtained from Figure
86, is converted to intrinsic permeability by the relationship
C« K v/g (93)
where g = gravitational acceleration. Thus:
c « (0.107 cm/min) l19^^)"7^/^] = 2.04 x io~8 cm2
Substitution of the above values into equation (90) yields
 n 12 x lo'ifflinf1'0 x 10"3C2.Q4 x loyn
 (1.12 x 10 )(88)^ i.i2 x io« J
= 1.98 x 10"5 cmVsec = 1.19 x io"3 cmVmin
The desired coefficients AC, BC, CC and DC are computed from equations
(77), (78), (79) and (80), respectively, as follows:
ACio _ (M6 } + [1.19 x 1Q3] _
A(\  [2(10)J I io2 J " °*0030119
216

0.25  0.20
I 10 I
0.25
10
cc10 =
tf
.06 ] fl.19 x lp3<
(10)J + I 10*j
= 0.0030119
DC
.. = o * 0 * [t°25ff'25?) . 0.13
The finitedifference equation for j = 4 and n = 10 is, therefore
 (0.0030119)c10 + (0.0750238)c10  (0.0030119)c10 = (0.13)
The system of equations involving all nodes is solved in a manner
analogous to that in the flow model.
PROGRAM DESCRIPTION
Quantity Portion
Summary—
The purpose of the model is to solve numerically the onedimensional,
unsteadystate flow equation for an unsaturated porous medium. Ini
tially, the program considers a given vertical column of soil with a
certain pressure (moisture content) profile. This column is parti
tioned into a discrete number of nodes, as shown in Figure 78. At
each time step, the program determines the values of hydraulic pressure
resulting from the values of two boundary conditions, namely: (1)
the combined effect of rainfall and evaporation at the top of the pro
file, and (2) net recharge rate at the bottom. Since these boundary
conditions change with time, it is necessary to conduct an inter
facing with the overland flow and groundwater programs.
In addition to the initial tension profile and boundary conditions,
some relationships involving moisture, tension pressure and hydraulic
conductivity must be provided as inputs to the model. These relation
ships are essentially fundamental properties of each soil. Therefore,
the program computations must be conducted independently for each major
217

soil type. Special formulas are also included to compute withdrawals
(sinks) of water at each node in the soil profile.
The output from the unsaturated zone flow model consists basically of
the values of moisture content and velocity, obtained from the corres
ponding pressure heads at the various nodes. One table is printed out
for every time step. The velocities affect the movement and disper
sion of substances dissolved in the water.
The computer model consists of a main program (designated as MAIN)
and several subroutines: COEF, EXTRA, VFLOW, WTABLE, PSIVAL, PSIMID,
SINK and MATRIX. Also, four function subprograms are utilized to per
form moisturetensionconductivity conversions. The names of these
functions are: CSL, CK, TENS and HUMID. A summary of the main pro
gram, subroutines and functions follows.
Main Program—
The main program first reads in the data. A partial list of data is
given below:
1. Number of nodes and time steps.
2. Rates of rainfall, evaporation and net recharge flow at each time
step.
3. Vertical distance between node and interval between time steps.
4. Initial values of moisture content for all nodes.
Subsequently, some of the values above are converted into consistent
units. In the first time step, the moisture values given are used
to evaluate the various soil moisture functions.
In the second time step (N =2), the information from the previous step
and the boundary conditions are utilized to compute the values of
pressure at each node. In carrying out these computations, MAIN calls
subroutines PSIMID, COEF, PSIVAL (in that order) and various function
subprograms. Thereafter, it calls subroutine VFLOW to determine
the vertical velocities at each node in this time step.
For the following time steps (N > 2), the calculations are slightly more
involved than those for N = 2. Subroutines EXTRA and WTABLE are called
in addition to those above. Tables of moisture values for each node are
printed.
At any time step, the option of plotting the moisture versus depth
profile can be exercised. This can be accomplished by calling a spe
cial subroutine package.
Subroutine COEF
At the jth node in the soil profile, it is possible to write an equa
tion of the form of (55) which relates the pressure head, \p, at the
218

node in question to the t»'s at its adjoining nodes. The coefficients
associated with these three heads are designated in the program as
A?, B? and C?, in which the superscript n denotes the time step. The
coefficient designated as D? contains all known constants in the equa
tion. In a soil profile containing L nodes, it follows that L equa
tions of the general form described above can be written. The solu
tion of the resulting system of equations (22) requires that L values of
the coefficients A., B., C, and D1? be calculated at every time step.
Subroutine COEF performs the calculations for these coefficients.
Two nodes, one lying on the upper boundary and another on the lower
boundary of the soil profile, have only one adjacent node instead of
the usual two. The calculation of the coefficients associated with
these two boundary equations is therefore a special case. It consists
of the application of Darcy's law between the boundary node in question
and its adjacent node. The flow between the nodes is set equal to the
known boundary condition flow.
The values of A1?, B., C. and D1} at each node obtained by subroutine
J «* J J
COEF are used elsewhere in the program to determine the tension head at
every node.
Subroutine EXTRA
To compute the coefficients mentioned above in subroutine COEF, it is
necessary to determine first the value of some moisture properties of
the soil, such as hydraulic conductivity, at every node. Unfortunately,
these properties are a function of the respective pressure heads, which
are not yet known for the time step in question.
For this reason, these pressure heads most be predicted using their
computed values from previous time steps. Subroutine EXTRA performs
these extrapolations. It may be noted that EXTRA is not called at
219

time step N « 2, since only one set of previous pressure heads is
available at that time.
Subroutine PSIMID
To determine the coefficients An, B1?, C? and D1? in subroutine COEF,
it is first necessary to evaluate certain soil moisture functions, such
as conductivity. The numerical technique requires that some of these
functions be evaluated, not at the nodes themselves, but rather at
points midway between the nodes. Subroutine PSIMID computes the
values of pressure head, r), at such midpoint locations, using linear
interpolation.
Subroutine SINK
The option exists in this program to specify certain rates of with
drawal of water (sinks) at each node in the soil profile. These with
drawals can consist of water uptake by plants and/or artificial drain
age. In the former category, the uptake rates must be determined from
climatic and agricultural information. In the latter case, any one of
several semiempirical formulas to compute drainage rates may be used.
Additional information on these two types of calculations is given in
Appendix E.
Subroutine SINK adds the plant and drainage withdrawals for each node
in the soil column. The resulting values of total sink flows are stored
and then utilized in subroutine COEF to calculate the D" coefficients.
Subroutine PSIVAL
The purpose of subroutine PSIVAL is to assemble the information needed
to compute the values of pressure head, \p, at each node, and then to
call subroutine MATRIX to carry out the computations.
220

First, PSIVAL determines whether the top node is saturated or not.
Saturation at this node can be caused by two conditions: (1) net
rainfall rate (rainfall minus evaporation) greater than the saturated
conductivity of the soil, or (2) existence of a ponding situation due
to high previous rainfall. If either one of the above conditions holds,
the definition of the boundary condition at the top of the profile is
achieved by setting the pressure head at the top node equal to the pond
ing depth. If these conditions do not hold, the top node is not satu
rated and the usual upper boundary condition involving Darcy's law
applies.
After the proper boundary condition above has been selected, subroutine
MATRIX is called to compute the values of pressure head at each node.
Subroutine PSIVAL then predicts the depth of ponding (if any) for the
next time step.
Subroutine MATRIX
This subroutine assembles the previously determined values of A1?, B? and
C? into a lefthand side matrix, and the values of D1? into a righthand
side matrix, as shown in system of equations (22). It then proceeds
to solve for the unknown values of pressure head, ty» at each node.
As mentioned previously, the lefthand side matrix is of tridiagonal
form and therefore the system is most easily solvable by an elimination
procedure. The specific technique used herein is Gaussian elimination.
It may be noted that the occurrence of a ponding condition at the top
of the profile introduces some changes in the calculations.
Subroutine WTABLE
This subroutine determines the elevation of the water table at each time
step. To achieve this, it checks the previously computed values of
pressure head, starting at the bottom node and proceeding upward. After
221

it locates the first unsaturated node (pressure head less than zero),
it computes the water table elevation.
Subroutine VFLOW
Subroutine VFLOW first utilizes the values of pressure head, determined
by subroutine MATRIX, to find the total head (elevation plus pressure)
at each node. It then determines the flow of moisture at each node
from Darcy's law.
Plotting Subroutines Package—
This package is composed of subroutines GRAPH, PINE, CURVE, PPLOT and
function ACOS. The option to plot any tabulated results may be exer
cised from the main program by calling subroutine GRAPH. This subrou
tine, in turn, controls the flow of information between the other sub
routines in the package.
Function Subprograms—
Four function subprograms are used in this computer model. They all
involve basic soil moisture relationships, and their names and purposes
are:
1. Function CSL. Computes the slope of the moisture versus pressure
head curve, as a function of the latter.
2. Function CK. Finds the hydraulic conductivity for a given value of
pressure head.
3. Function TENS. Converts values of moisture content into pressure
heads.
4. Function HUMID. Converts pressure heads into moisture contents.
The above relationships may be determined in the field for any given
soil type, or may be approximated using literature data for similar
soils. The latter approach was used herein. Segmented straight lines
were used to fit nonlinear relationships.
222

Quality Portion
Summary—
The main purpose of SOILQUAL is to compute, at every time step, the
concentrations of any given substance at the various nodes in a soil
profile. The required input data for this program can be broken down
into two basic types. The first category consists of water quantity
and flow data, obtained by interfacing with the unsaturated flow pro
gram. These data include values of flow velocity and moisture content
at every node, in addition to initial and boundary conditions. The
second type is composed of water quality information, such as fertili
zer availability and rate coefficients for source and sink reactions in
the soil.
The two types of data discussed above are used by the program to compute
the values of coefficients ACn, BC1?, CC?, and DC? in equation (75) at
•^ J J J
every node. Some of these coefficients vanish in the case of the top
and bottom nodes, which have only one adjacent node instead of the
usual two. The program then proceeds to solve a system of simultan
eous equations to obtain the desired concentrations at every node.
The computer model consists of a main program, designated as MAIN, and
subroutines GROUP, MAT2 and OUTRIT. Following is a discussion of
both the main program and the subroutines.
Main Program—
Initially, MAIN sets up data storage blocks from which to retrieve
water quantity information previously obtained by the flow program.
Other blocks are designated to store water quality data read into the
model, and other water quality parameters generated by the model itself.
223

After reading in the data, the program performs some simple calcula
tions involving conversion factors, to obtain the desired consistency
in units. Next, it utilizes hydraulic results from the flow program
to compute the value of the dispersion coefficient at every node.
Subsequently, a loop involving all time steps carries out the follow
ing instructions, in order: (1) calls subroutine GROUP to assemble at
every node the various terms in equation (74) into an equation of the
form of (75), and (2) calls subroutine MAT2 to perform the matrix
operations required to solve for the concentrations of substances at
each node. At the conclusion of this time loop, MAIN calls subroutine
OUTRIT to print out all the pertinent results.
Subroutine GROUP
At every time step, GROUP determines the values of coefficients AC, BC,
CC and DC for every node in the soil system. Recall that the top and
bottom nodes constitute special cases. The structure of this subroutine
is very similar to that of subroutine COEF in program INFILT.
Subroutine MAT2
This subroutine inserts the values of AC1?, BC1? and CC1? made available
by GROUP into a lefthand side matrix. Also, it inserts the corres
ponding values of DC. into a righthand side vector. Subsequently it
solves the system of equations in a stepwise fashion using Gaussian"
elimination. It may be noted that subroutine MAT2 is a modified
version of subroutine MATRIX in the unsaturated zone flow model.
Subroutine OUTRIT—
After subroutines GROUP and MAT2 have been called for all the time steps
in the period of study, MAIN calls subroutine OUTRIT to print out the
results. The most significant information provided by OUTRIT is a table
of concentration values at all nodes for all time steps.
224

APPENDIX C
DETAILS OF GROUNDWATER MODEL
INTRODUCTION
Shown herein are sample calculations for the groundwater flow model
as described in Section VI. Space limitations preclude the presenta
tion of counterparts for the quality formulation (Bedient, 1972).
Also included are descriptions of the subroutines in both the f;ow
and quality models. Consistent English units are used. Inasmuch as
the data are hypothetical, equivalent metric units are not given.
FLOW CALCULATIONS
The basic finitedifference relationship (Freeze, 1966) used in the
calculations expresses the hydraulic head at a given node in terms of
the heads at its adjoining nodes and other constants:
aKvCl,J)cf)(I,J+l)]
+ ccKv(I,Jl) + aK
where = hydraulic head, ft
K, , K = horizontal and vertical permeabilities, gal per
* v day/ft2
o
 (dimensionless)
Ax * horizontal spacing, ft
Az = vertical spacing, ft
225

Any consistent permeability units could be utilized in this formula,
since Kh and Kv act as weighting factors.
To illustrate the usage of equation (94), it is possible to isolate a
portion of a grid system, as shown in Figure 81 along with assumed
values of the permeabilities and geometric spacings. Also shown are
the heads at the various nodes. A new value of head at node (I,J) can
be determined by substitution into (94), as follows
(20) (70.0) + (20)(74.0)
+
'100'
50
V J
(40) (69.1) +
20 + 20 +
100
SO I
« f
rioo'
L50J
(40) (75. 0)1
(40) +
100
[50
^ /
(40)
= 72.2 ft
Next, a relaxation factor is applied to the computed value to promote
convergence in the calculations. The formula used for this purpose is
= wn(I,J)
(I.J)
(95)
in which w denotes the relaxation factor. The subscripts p, n and c
indicate the previous, new and corrected values, respectively. Expres
sion (95) above is the same as equation (39). In this study, the value
selected for w was 1.85 (Freeze, 1971). Substitution of the appropriate
values into equation (95) yields
(I,J)  (1.85) (72.2) + (0.85) (71.1) = 73.2 ft
Equations (84) and (95) are applied to the various nodes in the
groundwater system in an iterative fashion. It may be assumed that
after all iterations have been completed, the values of head (ft)
at the five nodes in question are:
226

COORDINATE AXES:
IZ
Ni
fO
HORIZONTAL AND VERTICAL
PERMEABILITES, gol/doy/sq ft
FOR ALL NODES:
Kv=40
VALUES OF HEAD ft AT
nth ITERATION'.
(II,J) = 70.0
4>(I,J) = 71.1
0(I + I,J) = 74.0
(I,Jl)= 69.1
0(I,J*I)= 75.0
Figure 81. Isolated portion of a groundwater grid system, showing
geometry and assumed heads and permeabilities.

(Il,J) • 71.5
(I,J) * 73.0
<>(I+1,J)  74.4

and PREDIS. Following is a description of both the main program and
its subroutines.
Main Program—
The main program first reads in the data, some of which are listed
below:
1. A matrix (NCLASS) for node classification. The numbering system
indicates, for instance, whether the node lies on the water table or
another boundary, or whether its head value is known. This matrix
remains unchanged throughout the program,
2. Horizontal and vertical permeabilities for every node,
3. The values of head at nodes where it is known, and trial values
at other nodes.
The program then proceeds to write out the matrix PHI, which contains
the information in item (3) above. Next, the iterative procedure to
calculate the head at each node where it is unknown is begun. The head
at a given node is expressed in terms of the heads and permeabilities
of its surrounding nodes, its own permeability and the horizontal and
vertical spacings. However, the exact form of the expression varies
with the node type as designated in NCLASS. For this reason, several
conditional statements in the program check the value of NCLASS for the
node under consideration to ascertain the proper formula to be used.
Most nodes are of the interior type, that is, surrounded by four ad
jacent nodes, and the standard formula for computing head is applica
ble, Other nodes lie on noflow boundaries of the grid system and are
thus surrounded by fewer than four nodes. In this case, subroutine
BOUND is called to perform the! head computation. If a certain discon
tinuity effect occurs at the water table, the head at the affected
node is determined by subroutine TABLE. In all the above cases, a
relaxation factor is applied to the value of head following its cal
culation.
After the iterative process to determine the heads throughout the grid
system has been completed, MAIN prints out these values in tabular
form. It may be noted that any grid nodes located above the water
table lie outside the groundwater system. For identification purposes,
their head values are printed out as zeros.
Next, subroutine EQUIP is called to determine the location of points
of equal head. Subsequently, MAIN calls subroutine VEL to determine
the velocities throughout the grid system. The last subroutine called
is PREDIS. Its function is to prepare the results obtained for use as
input to the water quality program.
229

Subroutine BOUND—
The subroutine was designed to compute the head at nodes located on a
noflow boundary. The existence of such a boundary can be brought
about by an impervious formation or by the flow patterns of a basin.
From Darcy's law, it follows that the hydraulic head on both sides of
a noflow boundary must be equal.
Since nodes lying on boundaries are surrounded by fewer than the usual
four nodes, BOUND must first establish the location of the missing
node(s) (up or down, right or left) with respect to the boundary node
in question. Subsequently, it implicitly creates an "image" or imagin
ary node at that location. To satisfy the noflow condition, the head
at this image node is set equal to that of the node located directly
opposite it across the boundary. Finally, the computation for the head
at the boundary node in question is carried out,
Subroutine TABLE—
In the unlikely event that a node lying to the right or left of an
interior node is located above the water table, the head computation
becomes another special case. Subroutine TABLE first establishes which
of the adjacent nodes (left, right, or both) lie above the water table.
It then branches out to perform the appropriate calculation.
The general condition described above is brought about by large varia
tions in the slope of the water table. By decreasing the spacings in
the grid system, the need to utilize subroutine TABLE can be minimized.
Subroutine EQUIP
This subroutine determines the location in the grid system of points
having certain values of head. It uses linear interpolation to achieve.
this purpose. A table of results is printed.
Subroutine VEL—
Subroutine VEL uses Darcy's law to compute the horizontal and vertical
velocity components throughout the grid system. It then uses these
components to determine the resultant velocity vector at each node.
The direction of this vector is computed as an angle, in degrees,
measured counterclockwise from the positive horizontal axis. Both
the magnitude and direction of the velocity vectors are printed out
in tabular form.
Subroutine PREDIS
The function of PREDIS is to assemble the results obtained throughout
program PHI for subsequent transfer to the groundwater quality node.
The variables used in the interfacing are the velocity components,
horizontal and vertical spacings, time step and others. This inter
facing can be carried out using punched cards or by means of common
storage locations.
230

Quality Portion
The computer program used herein has basically the same sequence of
operations as that developed by Shamir and Harleman (1966), with some
modifications introduced by Bedient (1972). Since the latter has
provided a detailed description of the program, only a brief summary is
given below.
The computer model consists of a main program (MAIN) and the following
subroutines: CINPUT, TRIX, TRIY, EXPLIC, COEFX, COEFY, CFEXP, VCALC
and OUTPUT. The function of MAIN is to read in preliminary data,
compute the dispersion coefficients, and control the sequence of opera
tions performed by the various subroutines.
Subroutine CINPUT reads in concentration data used as initial and
boundary conditions. The purpose of COEFX and COEFY is to provide the
values of the coefficients and righthand side constants for all equa
tions of the type of (40) and (41), respectively. The two systems of
equations corresponding to these types are solved by TRIX and TRIY.
After the computations involving the alternating direction implicit
procedure have been completed, subroutine EXPLIC obtains the final
concentration values using an explicit scheme. The coefficients used
in this matrix solution are provided by CFEXP.
If desirable, it is possible to read in velocity values from a tape
by calling VCALC. This subroutine also performs some necessary aver
aging calculations. Finally, subroutine OUTPUT prints out a table of
pollutant concentrations at each grid node, at any given time step.
231

APPENDIX D
DETAILS ON INTERFACING METHODOLOGY
A generalized schematic representation of the overlay procedure used
in this study is provided in Figure 82 which shows the basic branching
system containing the main program and various model subroutines. Each
of these model.ysubroutines corresponds to a major partition in the
hydrologic system and, in turn, contains other subroutines of lesser
hierarchy.
In the present formulation, any given subroutine can only be called
once during a simulation run. For instance, subroutines SUB1A and SUB2A
might in effect be identical, but their names cannot be the same. The
capability of reusing a subroutine name at different times in the
simulation would have required additional job control language (JCL)
statements and input/output instructions. These complexities were
thought to be unwarranted from a demonstration standpoint. Also, the
current relative inefficiency of the computations can be reduced by
the use of object decks, preferably within a program library.
Example JCL instructions for the overlay procedure applicable on the
IBM 370/165 at the University of Florida are given in Figure 83. The
overlay requires that the program setup be carried out through the
linkage editor, instead of the less sophisticated loader. Source or
object programs may be used, but their locations in the deck differ.
The different data sets are defined toward the end of the deck. In
the University of Florida system, data sets 5, 6 and 7 are reserved
for reading, writing and punching cards, respectively.
Lists of interfacing variables between the different types of model
subroutines are provided in Figure 84. The pertinent information for
the Stormwater Management Runoff and Receiving Models is available
elsewhere (Environmental Protection Agency, 1971).
A listing of the main program used for an overlay run may be found
toward the beginning of Appendix G. Defined therein are the names of
model subroutines and data sets used for the various subsystems of the
Lake Apopka basin. Also included is the calling sequence for the
model subroutines.
232

MAIN
PROGRAM
to
SUB I
SU8IA
SUB IB
SUB 1C
SUB 2^
SUB2A
SUB2B
SUB2C
MODEL SUBROUTINE
SECONDARY SUBROUTINE
SUB 13
SUB I3A
SUB I3B
SUB I3C
COMPUTER TIME
Figure 82. Generalized hierarchial representation of overlay procedure.

ro
04
// EXEC F4GCLM , PARM.LKED = 'LIST, MAP, VLY'
// F$RT. SYSIN DD *•
(Insert Source Programs)
/*
// LKED. SYSIN DD *
(Insert Object Programs) Arbitrary Name
OVERLAY ALP
INSERT SUB I, SUB IA , SUB IB , SUB 1C , ...
OVERLAY ALPHA
INSERT SUB 2 , SUB 2A , SUB 2B , SUB 2C,
OVERLAY ALPHA
INSERT SUB I3,SUBI3A,SUBI3B,SUBI3C,
/ K sData Set Number
/*"•*"
GO. FT01F001 DD UNIT = SYSDA, SPACE = (CYL ,{ 2,1))
GO.FT02F001 DD UNIT = SYSDA, SPACE = (CYL , (2 ,1))
//GO. SYSIN DD*
Figure 83. Job control language (JCL) instructions for
overlay procedure valid on IBM 370/65 computer
at the University of Florida,

DATA SET
From Model
Unsaturoted
Zone
Time Step
No. of Time Steps
Rainfall Rates
Infiltration Rates
Ponding Depths
To Modeh
Kinematic Wave
Overland Flow
Kinematic Wave*
Overland Flow
Time Step
No. of Time Steps
Flows
Nutrient Fluxes
Receiving
(Lake Apopko)
Stormwoter Management
Runoff
Time Step
No, of Time Steps
FJows in Mole Drains
Nutrient Fluxes through
Mole Drains
Receiving
(Muck farms canals)
Unsaturated
Zone
Water Table Elevations
(heads)
Nitrate Concentrations
Groundwoter
Figure 84, List of major interfacing variables for the
various models.
235

APPENDIX E
DETAILS ON APPLICATION TO THE
LAKE APOPKA BASIN
INTRODUCTION
This appendix contains brief discussions on hydrologic and nutrient
budgets conducted within the Lake Apopka basin. These are followed
by explanations of various agricultural calculations. Finally, several
lists of data used in this study are provided.
HYDROLOGIC AND NUTRIENT BUDGETS
Budgeting efforts have been attempted for two portions of the study
area: (1) Lake Apopka and (2) the muck farms. Calculations for the
former were carried out as a supplement to the more detailed modeling
computations which constituted the main thrust of this study. The
latter budget was conducted by Heaney, Perez and Fox (1972) for the
East Central Florida Regional Planning Council.
Lake Apopka Budget
Numerous sources and sinks of water and associated nutrients were
considered in the analysis. For some of the budget items, such as
Gourd Neck Spring and the ApopkaBeauclair Canal, the available data
were adequate. Most of the data deficiencies involved items related
to the muck farms, namely: (1) pumpage rates, (2) flows through the
flap values, and (3) seepage through the earthen dike.
The approach taken in the hydrologic budget was to lump the three items
above into an unaccountedfor inflow term. Monthly estimates of other
sources and sinks, along with lake stage data, were used to solve for
the unaccountedfor inflow. The results obtained for years 1961 through
1970 indicated the existence of significant average annual flows from
the muck farms to the lake. However, in some months the computed
flows occurred from the lake to the farms.
Measured nitrogen and phosphorous concentrations at the various sources
and sinks were utilized in the estimation of a nutrient budget for the
236

lake. Due to lack of data, nutrient cycling between the lake sediments
and the water was not considered in the calculations. Also, the Hydro
logic indeterminacy involving the three items at the lakemuck farm
interface posed severe difficulties in computing a net nutrient flux.
At best, only limiting values could be estimated. In any case, the
various budget results indicated a significant overall net influx of
nutrients into Lake Apopka.
Space limitations preclude the presentation of the specific assumptions
and results of both the hydrologic and nutrient budgets. However, this
information has been tabulated and filed for possible subsequent use.
An examination of the hydrologic budget developed by Anderson (1971)
should also be useful.
Muck Farms Budget
The hydrologic indeterminacy mentioned above, among other data defi
ciencies, precluded an accurate estimate of both the hydrologic and
nutrient budgets within the muck farms. Since the overall results
may be found in a report by Heaney, Perez and Fox (1972), it is more
pertinent to discuss certain items which provided insight into modeling
approaches.
One of these items consisted of the nutrient contribution by rainfall.
It was determined that this source is significantly smaller than ferti
lizer application. Therefore, rainfall was not considered in defining
the upper boundary conditions in the quality portion of the unsaturated
zone model. This omission was also made in the sandy soils area.
Computed yearly average fertilizer application rates were used to
estimate nutrient concentrations at the uppermost node of the soil
profile. Estimated harvest weights and nutrient contents of the var
ious vegetable crops in the muck farms were utilized to determine nu
trient withdrawal (sink) rates from the soil.
Finally, consideration was given to the release of nutrients by the
muck soil. If the soil losses in the muck farms were comparable to
those in the Florida Everglades of 1 in/yr (2.5 cm/yr) (Stephens,
1969), the corresponding nutrient losses could be very significant.
CALCULATIONS
Agricultural calculations conducted in the Lake Apopka basin involved
both types of overland flow models and the unsaturated zone formulation.
237

Overland Flow
The usage of the duBoys formula (Einstein, 1964) to compute sediment
loads at the catchment outlets was illustrated in Appendix A. It is
also of interest to consider the routing of phosphorous and nitrogen
in overland flow.
Certain soils (e.g., muck) contain large amounts of nutrients, pre
dominantly in organic form. Thus, these nutrients are transported
as part of suspended solids loads. Due to the short duration of over
land flow phenomena, it was considered unlikely that significant amounts
of nutrients would be released by the soil particles. However, it is
expected that significant releases could occur in the receiving bodies
of water.
Since organic N and P fluxes can be estimated from soil chemical data
and predicted suspended soilds loads, separate attention was given to
the inorganic forms of these nutrients.
Phosphorous Routing—
The only inorganic phosphorous source considered herein was fertilizer
application. Data on fertilizer loading rates for the study area were
coupled with generalizations on phosphorous behavior in soils (Stanford,
England and Taylor, 1970).
Agricultural phosphorous is fixed rapidly to most soils, and thus tra
vels with suspended particles in overland flow. For this reason, the
weight fraction of this substance in the soil was sought.
The calculations may be illustrated by referring to Figure 85. Part
(a) of this figure constitutes a top view of an arbitrary square plot
within the sandy soils region. The area of this plot is one acre.
Part (b) is a lateral crosssection of the plot. Illustrated is the
estimated phosphorous loading rate corresponding to fertilizer appli
cation in citrus groves of 60 Ib/ac/yr (64.4 kg/ha/yr), obtained from
the East Central Florida Regional Planning Council. Since the citrus
groves occupy approximately onehalf of the total area of sandy soils,
the areal average application rate was assumed to be 30 Ib P/ac/yr
(32.2 kg/ha/yr).
Several simplifying assumptions were necessary to estimate the distri
bution of the applied P within the soil profile. One was that 90
percent of the applied amount remained within the top six inches of the
profile. It was also assumed that the distribution was uniform with
depth within this range (see Figure 85).
Following application, phosphorous is slowly removed from the soil by
plants. The remaining amount at any time depends, therefore, on the
238

K>
04
AREA = 1 acre
r
•O
'
"1
CITRUS GROVES
PART A PLAN VIEW
P APPLICATION RATE
3O
M ll?l
G1
P CONCENTRATION
. PROFILE I
WATER TABLE '
v
PART B LATERAL VIEW
Figure 85. Illustration of fertilizer application on a plot in the
sandy soils.

time elapsed since application and the growth characteristics of the
crop. A gross assumption was made that, on a timeaverage basis, the
P remaining in the soil was onehalf of the applied amount.
The concentration of P within the top six inches (15 cm) can be com
puted as:
(0.90) [IJ(30 Ib)
°P ~ f i f+\ ( f**}
= 6.18 x lo'" = 9.93 x
About 43 percent of the sandy soil is occupied by voids. Thus, the
concentration in the soil itself becomes
'  100 r  1 f)ft x in3 Ib P _ i 77 x ln2 Kg P
Cn TTnTrTTpD ~ 1'U° X JU •——= = 1.73 x 10 —s—*
°P (10043) P ft3 sand m3 sand
That is, suspended solids fluxes must be multiplied by c to obtain
the phosphorous fluxes.
The calculations for the muck soils were based on similar assumptions.
The value of Cp obtained was 4.25 x 10"3(lb P/ft3 muck); or metric 
6.92 x 10'2(Kg P/m3 muck).
Nitrogen Routing—
According to Stanford, England and Taylor (1970), nitrogen applied in
fertilizers is generally readily converted to nitrate ion. This sub
stance leaches rapidly through soils.
In this study, the contribution of agricultural nitrogen to overland
flow was viewed to occur, not as part of an erosion process, but rather
as diffusion upward from the soil water. It was thought that the rate
of transfer might decrease exponentially ^ith time, and such a relation
ship was incorporated into the model. However, due to lack of data
these transfers were tentatively assumed to be zero.
240

Unsaturated Zone
Calculations for the unsaturated zone involved consideration of diverse
items, ranging from soil curves to source/sink effects within the soil.
Hydraulic Conductivity Curves—
Plots of assumed pressure versus conductivity relationships for sandy
and muck soils are provided in Figure 86 and 87, respectively. Only
the saturated conductivities for these two soils could be estimated
with reasonable accuracy (Stewart, Powell and Hammond, 1963; Harrison
and Weaver, 1958). No data were.available for unsaturated conditions.
In the muck, this deficiency is less likely to be severe because the
soil is maintained close to field capacity. The curve for sandy soil
was set by visual comparison to curves reported by Freeze (1969).
Mole Drains—
The expression used to compute the flows through the mole drains in
the muck farms was the following (Harrison, 1960):
4ifJ2d* + hl
Q . 15iL L (96)
where Q = discharge, m3/hr per unit length in m
K = hydraulic conductivity, m/hr
h = vertical distance between water table and drain, m
do  thickness of an "equivalent layer", m
e
S = spacing, m
The discharge on a per unit area basis, q, can be obtained by dividing
the above expression by the spacing, as illustrated in Figure 88,
part (a):
.9.
S C97)
Accordingly, the units of q are m3/hr/m2, or m/hr. For consistency
with equation (50) in Appendix B, the units of q are converted to cm/mini
241

10
ISO
E
u
100
tu
o:
en

K>
30O
o
^200
tr,
^
CO
w
LU
£100
0.01
MUCK SOILS
S»vsK
0.1 0.2 0.3 03275
HYDRAULIC CONDUCTIVITY, K, cm/min
0.4
Figure 87,
Assumed pressure versus hydraulic conductivity
relationship for muck soils.

NJ
•MOLE DRAINS
PART A PLAN VIEW OF UNIT ELEMENT
rj[
I ',!
X X X
CAN At. A
II
I
t_JUNCTION 
AN
ll
II
I
1,1
1
!!
OUTLINE OF CATCHMENT
PART B TOP VIEW OF THE MOLE DRAINS
AND CANAL SYSTEMS
Figure 88. Illustration of layout of mole drains.

100 ^
6o2hT
Finally, the discharge must be expressed on a per unit depth basis:
in which Az represents the vertical spacing (cm) . The variable q .
(1/min) represents the sink term for mole drains. d
To compute the contribution of the mole drains to a given junction in
the receiving model, it is necessary to determine the total length
of drains feeding into such a junction. An approximate relation for
this length, L, is:
(99)
in which A represents the area of the runoff catchment and S denotes
the spacing. The total flow then becomes
QL (ioo)
where Q is the flow per unit length computed from equation (96) •
Appropriate unit conversions must be carried out in equation (100).
The nutrient fluxes into the receiving junctions were computed using
the flows from equation (100) and concentration data collected by
Hortenstine and Forbes (1971) .
Top Boundary Condition—
The concentration value at the top node, resulting from fertilizer
loads, depends upon the time elapsed since application, uptake by plants
and other conditions in the soil. For simplicity, the computation of
this value is illustrated for the situation immediately following
application.
245

The loading rate for phosphorous in the sandy soils was estimated at
30 Ib/ac/yr (32.2 Kg/ha/yr) . It is assumed that following application,
the entire amount becomes uniformly distributed with depth, down to
a distance Az/2 from the ground surface. The vertical spacing, pre
viously selected, was 17 cm. For convenience, a saturation condition
may be assumed to exist at the top node. Accordingly, the volumetric
moisture fraction is 0.43.
The concentration in the soil water, c, can be computed as follows:
30 Ib x 454^ x I03m&
_ weight m _ IP _ g _
volume  0>43 x 1 ac x 43<560 £_ x 8>5 cm x 1 J£ x 2g>32 i
elC • OU Cm XL
= 9.2^
Analogous calculations may be conducted for nitrogen.
Plant Withdrawal s
Quantity— Data on evapotranspiration rates for a given area provide
the starting point for calculations involving water uptake by plants.
For instance, for the citrus groves in the Lake Apopka basin it was
estimated that four fifths of the total water loss from the soil con
stituted transpiration by the crops. Data on total evapotranspiration
rates for citrus groves in Florida have been compiled by Hashemi and
Gerber (1967).
It is known that water uptake by crops is not uniform throughout the
root zone but, rather, decreases with depth. Several distributions of
this uptake were considered in preliminary calculations. However, it
was found that overall model results for short simulation periods were
insensitive to plant withdrawals, which are significant only on a long
term basis. For this reason, it was decided to assume a uniform with
drawal pattern and thus simplify programming instructions.
To illustrate the computation of the sink terms, it may be assumed that
the total uptake by citrus is 0.10 in/ day (0.25 cm/ day) . Assuming also
a root zone depth of five feet (1.5 m) and the designation of 10 nodes
within that zone (see Figure 89) the withdrawal rate at a given node is '
W = 0'?° in/day = 102 in/day = 1.76 x 105 cm/min
ID nodes
246

K)
JUPTAKE RATE = 0.10in /day
Figure 89. Illustration of water uptake by citrus.

This value must be converted to a per unit depth basis:
1.76 x io5 SL
w = ~ =  SH>  = 1.03 x io6 4
Az e t+ v in A cm v 1 mm
5 ft x 30.4 x
The units of w are consistent with those specified for a sink in equa
tion (50).
In the computer program, the soil tension at each node is checked
prior to computing the sink term, if the tension exceeds that for the
wilting point, a withdrawal of zero is specified. It is perhaps likely
that, in practice, a nowithdrawal condition might result from even
smaller tensions.
Quality— The nutrient withdrawal pattern within the soil profile should
follow that of water uptake, since nutrients are dissolved in water.
However, for simplicity, nutrient uptake was assumed to be uniform with
depth. Also, seasonal effects and crop growth patterns were not
considered.
In the muck farms, the nutrient uptake rates were estimated from the N
and P contents of the crops and their harvest weights. Yearly with
drawals were 96.0 Ib/ac (103 Kg/ha) for N and 8.9 Ib/ac (9.6 Kg/ha)
for P (Heaney, Perez and Fox, 1972) .
Calculation of the nitrogen sink term may be illustrated by converting
the withdrawal rate into equivalent units of nitrate
96.0 Ib N = 96.0 x = 427 Ib NOJ * 194 Kg NOJ
Assuming a root depth of one foot (0.3 m),and a node spacing of 6.7 cm,
the uptake for a given node becomes, on a per month basis:
u= 427 x x    ,or
'acyr 12 moj 1 ft x 30.5 cm >
U7.8
• acmo hamo
248

For consistency within equation (73), the desired units for this
withdrawal rate are (rag/ 1/min):
u =
/ \
7.8 ^
1 ra°J
mg
1 *\ 99 f* JllJtll
10.35 rr
ID
mo .
(1 ac)(6.7 cm) 4.05 x 10" —^—16
^ accraj
2.97 x IQ" mg/1
9 min
in which 6 is the volumetric moisture fraction. Assuming 6 = 0.75,
the uptake becomes
u = 3.97 x 10* mg NOT/1/min
The calculations for phosphorous follow the same steps.
Little information was available on nutrient uptake by the citrus in
the Lake Apopka basin. One approach used to estimate this uptake was
to examine the ratios of plant nutrient withdrawals to fertilizer load
ing rates in the muck farms and attempt a correlation. Data on this
aspect of citrus would be desirable.
Soil Source/Sink Effects
Since these soils are unlikely to release nutrients, the soil medium
should act as a net sink for both nitrate and phosphorous, especially
the latter. However, due to lack of data, these sink terms were tenta
tively set at zero.
On the other hand, muck soils are known to release nutrients. Little
quantitative information is available on the release mechanism, but for
illustrative purposes it may be tied to the consolidation of these soils
(Stephens, 1969).
The nutrient release resulting from a hypothetical one inch (2.54 cm)
yearly soil loss can be estimated from data on nutrient contents of
muck (Forbes, 1971). Since these data are expressed on a percent dry
weight basis, it is necessary to estimate the weight of one inch of
soil over an area of, say, one acre (0.4 ha): Assuming that the dry
muck unit weight is 100 lb/ft3 (1606 Kg/m3) the weight can be computed
as:
249

W = 1 inxL^lx 4,356 x 10" £ x 1 ac x 100 £r = 3.63 x lo^b
L£ 1T1 3C
 1.65 x 10s Kg
The nitrogen content of muck is 2.79 percent. Thus, the amount of this
substance released is:
r = x 3.63 x 105 = 10.1 x 103 = 10.9 x 103
100 ac ha
Assuming that this loss occurs over a period of one year and that all
the nitrogen is converted to the nitrate form
r  10.1 x 10' x = 4.48 x 10" „ 4.82 x 10*
acyr "tw* *w hayr
The units for this release rate are the same as the original units used
for the harvest removal estimates in the muck farms. Therefore, con
version into the desired sink units of mg/i/min is carried out in the
same manner as in the plant nutrient withdrawal calculations. Hope
fully, a better understanding of release mechanisms should suggest
improvements in the corresponding calculations.
LIST OF DATA
The scope of this investigation precluded a separate sampling program.
Nevertheless, substantial amounts of available data collected by various
agencies were assembled. These data are summarized in Tables 1 through
4. Indicated are the sampling frequency, agency and other pertinent
information.
Most of these data have been punched on computer cards. The original
tabulations, which are mostly in computer coding sheets, are also on
file.
250

TABLE 1
List of Hydrologlc Data for Lake Apopka Basin
fO
in
Item
Pumping rates
from muck, farms
to Lake Apopka
Rainfall
Stagearea data,
before and after
muck farms
Pan evaporation
Pan evaporation
Lake stage data
Rainfall data
Apopka— Beauclalr
Canal flows
Piezoicetrlc
levels
Location
Zellwood and
Duda muck farms
Clermont and
Lisbon stations
near Lake
Apopka
Lake Apopka
Lisbon station
near Lake Apopka
Lisbon station
Lake Apopka
Stations In
muck farms
North shore of
Lake Apopka
Orange County
and Lake County
Source
East Central
Florida Regional
Planning Council
(ECFRPC)
ESSA (Heather
Bureau)
uses
ESSA (Heather
Bureau)
ESSA (Weather
Bureau)
uses
Muck farms
uses
DSGS
Dates
19671970
19361970
19611970
19671970
19361970
19611970
19671970
19681971
Sampling Frequency
Dally, with gaps
Monthly and daily
. Monthly
Daily
Daily
Daily, with gaps
Daily
Twice per year
Comments
Hot broken down
by individual
pumps
Monthly 19361970
Daily 19671970
Areas available on
0.1 ft. stage
intervals

TABLE 1—Continued
to
tn
K>
Zten
Coordinates of
veils
Rainfall data
ApopkaBeauclair
Canal flows
Location
Orange County
and Lake County
ApopkaBcaucla.lt
Canal
North shore of
Lake Apopka
Source
uses
ECFRPC
uses
Dates
19361967
Sampling Frequency
Monthly
Monthly, with
numerous gaps
Consents
Longitude and
latitude

TABLE 2
U*t of Agrlcnltural Data for Lake Apopka Basin
Ite»
Land use data
Acres planted
with a given crop
Tertlllzer types
for various crops
Portions of faros
lying in each
runoff subcatchsent
Amounts of ferti
lizer applied to
•various crops
Areas of runoff
sub catchment
Farm areas
Location
Lake Apopka
has la
Muck ferns
Kuck faros
Muck, farms
Mock f aras
Mock farms
Mock farms
Source
Soil Conservation
Service (SCS)
East Central
Florida Regional
Planning Council
ECFKPC
ECFBPC
ECFEPC
ECFKPC
ECFSFC
Sates
Recent
19701971
19701971
19701971
19701971
19701971
19701971
Sampling Frequency
Weekly
COBEBditS

TABLE 3
Water Quality Data for the Lake Apopks Basin
(S3
tn
Item
Water quality
parameters,
Including N
and P
Water quality
parameters, as
above
Water quality
parameters, as
above
Water quality
parameters, as
above
Typical values of
various water
quality parameters
Various water
quality parameters
Ions and solids ,
Sediment core
samples including
Location
21 stations in
Lake Apopka
Shallow wells
around Lake
Apopka
Stations in
muck farms
canals
3 stations in
Winter Garden
Canal
10 stations in
Lake Apopka
21 sampling
stations in
Lake Apopka
21 bottom
sampling
stations in
Lake Apopka
Stations as
shown in
Source
Orange County
Pollution
Control Dept.
(OCPCD)
OCPCD
OCPCD
OCPCD
Florida State
Board of Health
(FSEH)
FSBH
OCPCD
EPA
(Schneider
Dates
19671970
1970
19681970
19691970
October 1962
and January
1963
First half
of 1964
Summer 1967
Report dated
1968
Sampling Frequency
Once every month
or 2 months
Monthly
Monthly
Every 2 weeks
Every 2 weeks
Once
Approximately 3
tines each station
Comments

TABU 3—Continued
Is)
VI
cn
Item
nutrients and
solids
Nutrient
concentrations
Effect of lake
bottom on fish
feed organisms
Biological
organisms
(mainly bottom)
Phosphorous data
Chemical and
biological
parameters
Nutrient
concentrations
Location
report
Rainwater,
shallow wells
and artesian
veils in Lake
Apopka basin
Lake Apopka
Lake Apopka
Muck farms
and Lake
Apopka
Lake Johns
(south of
Lake Apopka)
Percolating water
in muck faros
Source
and Little)
Appendix B
EPA
(Schneider
and Little)
Appendix C
FSBH report
dated 1964
OCPCD
Or. Forbes
(IFAS, Univ.
of Florida)
OCPCD
Drs. Forbes and
Hortenstlne
(IFAS, Univ. of
Florida)
Dates
Report dated
1968
19471963
19671970
19671970
19661969
1971
Sampling Frequency
Once
Every 2 years
Every 3 months
Monthly
Every 3 months
Seasonal
Comments
Qualitative
comments for
each station
Being prepared
for publication

TABLE 4
Miscellaneous Data Used in Computer Programs for Lake Apopka Basin
ro
tn
Item
Permeability data
Soil charac
teristics :
BOisture^tension
conductivlty
Geometry of
subcatchments and
canals . Station
locations. Link
age data.
Sub catchment
geometry and
linkage data
Fertilizer data
and ebil charac
teristics
Location
Vertical cross
sections in
Lake Apopka
basin
Muck soils and
sandy soils in
Lake Apopka
basin
Zelluood and
Duda muck
farms . Lake
Apopka.
Sandy soils
Muck and
sandy soils
Source
DSGS geologic
data
Various Soil
Conservation
Service data
Assorted naps
USGS quadrangle
maps
East Central
Florida Regional
Planning Council
Dates
Recent
Recent
Recent
Recent
Sampling Frequency
Comments
Used is ground
water model
Used in unsatu
rated zone
program
Used in runoff
and receiving
models
Used in runoff
model for undu
lating areas
(sandy soils)
Used in quality
portion of
unsaturated
zone program

APPENDIX F
SUMMARY OF ADDITIONAL REFERENCES
ON AGRICULTURAL NUTRIENTS
Although the state of the art in hydrologic modeling is at a relatively
advanced stage, much remains to be accomplished in the modeling of
agricultural water quality. Hopefully, the references mentioned herein
should be of some value in formulating needed improvements.
One source of nutrients not considered in this study was animal wastes,
due to their insignificance in the Lake Apbpka basin. However, this
source can often be important, as pointed out by Townshend, Reichert
and Nodwell (1969). These researchers provided figures on nutrient
contributions by different animals, expressed in terms of human popu
lation equivalents. They also discussed trends in animal populations
in Ontario, Canada. Webber and Lane (1969) studied the different
nitrogen compounds contained in manure and the ensuing reactions in
the soil.
Water quality characteristics of runoff from cattle feedlots were
investigated by Norton and Hansen (1969). Their study comprised both
a mathematical development of overland flow and a quality sampling pro
gram. Their results should provide useful insights into water quality
modeling approaches.
Other investigations on overland flow quality include those by Gilchrist
and Gillingham (1970) and by Nelson and Romkens (1970). These two
papers deal primarily with field studies on phosphorous movement re
sulting from different fertilizer application rates, and could be used
as guides for sampling programs in a given area. Data on the concen
tration of different ions in overland flow have been collected by Gburek
and Heald (1970). Bacterial levels in runoff have been reported by
Kunkle (1970).
The calculation of an inventory of nutrient sources and sinks in a
given area can be facilitated by the use of the tabular approach pro
posed by Schultz (1970). This was the basic approach used for the
257

budgeting calculations in the muck farms of the Lake Apopka basin.
Insight into the economics of fertilizer application is provided in
a paper by Fox and Allee (1970).
258

APPENDIX G
LISTING OF COMPUTER PROGRAMS
Included herein are listings for the following models: (1) unsaturated
zone, (2) kinematic wave overland flow and C3) groundwater. Listings
for the stormwatermanagement runoff and receiving programs are avail
able elsewhere (Environmental Protection Agency. 1971).
A listing of the executive program, used in connection with the overlay
technique, is provided first. Usage of this program is optional.
259

C
C
C
C
PAIN PROGRAM FOR INTERFACING
COKVCN BLOCKS ON?,. UNSAT1 AND NUM1 RENAMED ID TO SAVE CORE SPACE
CQMGN/DSETS/NGWIF,NUUIF1,NUNOV1,NUNRE1,NUNGWUNUNRE2»NUNIF2«
NU!40V2tNUN3V3iNUNOV

13
I*
15
16
17
18
19
20
21
NUNOV4
NUMF3
NUISGW2
NUNGW3
NOVRE1
NOVRE2
NRERE1
NRERE2
NOVRE3
UNZON3
UNZON3
UNZON2
UNZON3
OVLAN3
OVLArK
REC1
REC2
OVLAN1
OVLAN2
UN Z ON 3
GWBLOK
GWBLOK
REC1
REC2
REC3
REC3
REC3
c
c
c
c
c
c
c
c
c
c
END
SUBROUTINE UNZONl
COKMON/DSETS/NGWIP,NUNIFl,.NUNOVl.,NUNREl,NljNGWl,NUNRE;2,NUNIF2I,
1 MUNCV2,NUNOV3,NUNIOV4,NUNIF3,NUNGW2,NUNGW3,NOVREl,NOVRE2tNREREl,
2 NRERE2,NOVR£3
RETURN
ENC
SUBROUTINE UNZON2
CONMCN/DSETS/NGWIF,NUNIFl,.NLNOVl,NUNREl,NUNGWl,NUiMRE2,NUNIF2,
1 NUNOV2,NUNOV3tNUNOV4,NUNIF3,NUNGW2,NUNGW3,NOVREl»NOVRE2»NREftEl,
2 NRERE2.NOVRE3
RETURN
END
SUBROUTINE UNZON3
COMMON/OS ETS/NGW IF, NUNIFl,.NUNOVl,NUNREl,NUNGWl,NUNRE2tNUNlF2i
I NUNOV2tNUNOV3,NUNOV*»,NUNIF3,NUNGW2,NUNGW3NOVREl,NOVRE2tNREREl,'
2 NRERE2.NOVRE3
CALL UNFLO
CALL SCILQU
CALL IRAN
RETURN
ENO
SUBROUTINE UNFLO
C UNSATURATED FLOW (INFILTRATION) PROGRAM tNFL i
C CNECIMSNSIONAL (VERTICAL) UNSTEADY STATE INFL 2
C AT EACH TIME STEP, THE PROGRAM USES BOUNDARY CONDITIONS AND INFL 3
C INFORMATION FROM ThF PREVIOUS TIME STEPS TO COMPUTE THE TENSION INFL *
C HEAD, (PSD, AT EACH NODE IN THE VERTICAL PROFILE. INFL 5
C THEREAFTER, THE PROGRAM COMPUTES AND STORES MOISTURE FLUXES INFL 6
C BETWEEN NOCES^ HYSTERESIS EFFECTS ARE NOT CONSIDERED INFL 7
C .
C LIST OF VARIABLES
C NNODES'NO OF NOCES IN PROFILE
C NSTEPS»NO OF TIME STEPS
C TSTARTSTARTING TIME OF SIMULATION, IN MIN.
C TENO.ENDING TIME OF SIMULATION! IN MIN
C R.AINRAINFALL RATE
C EVAPEVAPORATION RATE
C GWFLUW«NET RECHARGE AT BOTTOM BOUNDARY
C DEL/.VERTICAL SPACING, IN CM.
C CELTTIME STEP, IN MIN.
C THENOTINITIAL MOISTURE CONTENT
C ZDRAINELEVATION OF ARTIFICIAL DRAIN
C ZBOTOMELEVATION OF BOTTOM NODE
C PSJNOTINITIAL PRESSURE HEAD VALUE
C PSI»PRESSURE HEAD
C CSLOPESLOPE OF MOISTURE VS TENSION CURVE
C COND*HYDRAULIC CONDUCTIVITY
C NOCSATNO OF UPPERMOST NODE THAT IS SATURATED
C A ACOEFFICIENT IN TRIUIAGONAL MATRIX USED TO FIND PSI VALUES
C B BCOEFF1CIENT IN TRICIAGONAL MATRIX USED TO FIND PSI VALUES
C C« CCOEFFICIENT IN TRIDIAGONAL MATRIX USED TO FIND PSI VALUES
C C« CCOEFFICIENT FOR RIGHT HAND SIDE OF MATRIX SCHEME USED TO FIND
C PSI VALUES
C PI,P2.P3,SYMBOLS FOR PSI AT MIDPOINT LOCATIONS (SEE FR£EZE,1969)
C C1,C2«CONDUCTIVITIES CORRESPONDING TO P1,P2
C REXCESS RAINFALL (RAINFALL MINUS EVAPORATION)
C CELPSIDIFFERENCE IN PSI VALUES
261

C ^AdV»ZQeL=eL6VAT!ON AT UPPER AND LOWER NODES, RESPECTIVELY
C PHIABV,PHIB2L=PRESSURE HEADS AT UPPER AND LOWER NODES,RESP£CTIVELY
C GRAC12=PRFSSURE HEAD GRADIENT BETWEEN NODES 1 AND 2
C ZTA3LE=ELEVATION OF WATER TABLE
C NNSAT*INCEX COUNTING MW3ER FOR NUMBER OF SATURATED NODES*
C STARTING AT BUTTOM
C HPCND=PONDING DEPTH AT GROUND SURFACE
C FSIGND=PSESSURE HcAD VALUE AT TOP NODE
C SATCCN=SATURATEC CONDUCTIVITY
C CPLANT=FLOW UPTAKE BY PLANTS AT A NODE
C CORAIN=FLOW WITHDRAWAL BY DRAIN AT A NOOe
c CSIC;K=TOTAL SINMWITHDRAWAD AT A NODE
C XL=ARRAY FOR CUEFFICIENTS A,B AND C (SEE ABOVE)
C RHS=RIGHTHAND SIDE VECTOR FOR 0COEFFICIENTS (SEE ABOVE)
C TEKS=TENSION PRESSURE
C CK=CONDUCTIVITY
C CSLSLOPE OF MOISTURE VS PRESSURE HEAD CURVE
C HUMID = MOISTURE CONTENT
C £TCRCP=WATER WITHDRAWAL(ET) BY CROPS IN STUDY AREA
C WILTK=WILTING POINT MOISTURE OF A GIVEN SOIL
C RINF=INFILTRATIQN RATE (DOWNWARD VELOCITY BETWEEN TWO TOP
C NODES) IN CENTIMSTERS/M1N
C ORGOTDEPTH OF ROOT £QN6 (FT)
C GWCLC=FLQW OUT CF A MOLE DRAIN (CFS/FT)
C JTA3=NOU£ AT WHICH WATER TABLE IS LOCATED
C HTAS=HEIGHT OF WATER TABLE ABOVE BOTTOM OF UNSATURATED PROFILEt
C IN F6ST,
COMMON/ 10 /PSINOT(50>,RAIN(20),EVAP(20) ,GV>FLOW(20) , ZTABLE120) i
1 PSI(5G,20)»M50,20)«B(50».2G),C(50»20) ,0150,20) ,E ( 50..20) ,F (50V 20)
2 CSLOPE(50,20), TH£TA( 50, 20) , THENOT1 50) , VELZ I 50, 20) ,PSI1 (50,20) ,.
3 PSI2(50*20>tPSI3(5C»20>f. NOOSAT(20),
4 CUNC(SO»20),PSIBOT(50)t"HPONO(501,QSINK(50i20)tQNET12(50>»
COMMCN/LAB/ TITLE ( 20 ) , XLABI 11 ) , YLAB ( 6) (HORI1(20) , VERT (6)
DIMENSION G1650)
DIMENSION TIT(20)
COM^CN/CROP/ETCROPtWILTM
COMMON/ROOT/ DROCT
COMKCN/MOLE/CMOLC150)
tINFL 9
INFL 10
IMFL 11
INFL 12
INFL 13
INFL I*
INFL 15
INFL 16
50)
COKMCN/PONO/DPONDf 50)
CQMMCN/T IKE/NT I fE,DTIME
COMKON/SOIL/KSOIL
COMMCN/LCCMOL/LfOLE
COM^ON/WTAB/ JTAB( 50 ) i HTAB( 50)
COMKCN/DSETS/NGWIF,NUNlFl,.NUNaVl,NUNREl,NUNGWlfNUNRe2,NlJNIF2.
1 NUrtOV2,NUNOV3,NUNOV4fNUNIF3fNUNGW2(NUNGW3,NOVREl»MOVR62iNRER
101
100
120
150
141
142
2 NRERE2.NOVRE3
READ (5tlOO) NNODES.NSTEPS
WRITE (.6,101) NNODES
FORMAT (' NO OF NQOes=«,I5)
FORMAT (213)
NTIME'NSTEPS
READ (5,120) TSTART.TEND
FORMAT (2F8.2)
RAINFALL IN INCHES/HOUR
READ 15,140) IRAIN(I),I»i,NSTEPS)
FORMAT (15F5.2)
EVAPORATION IN INCHES/HOUR
READ (5,150) (EVAP(I),I»1,NSTEPS)
FORMAT (15F5.2)
WRITE (6,141)
FORMAT C1',20X,'RAINFALL AND
WRITE (6,142)
lYT STEP''5X''RAIMFALL'*7X''EVAPORATION')
EVAPORATION IN INCHES/HOUR')
INFL 17
INFL 18
INFL 19
INFL 20
INFL 21
INFL 22
INFL 2»
INFL 26
INFL 27
INFL 28
INFL 29
INFL 30
INFL 31
INFL 32
*NFL H
INFL 3*
INFL 35
262

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
144
143
160
162
164
163
195
161
190
130
170
180
205
210
7500
7600
200
300
WRITE (6,144) I.RAIMI I),EVAP( I)
FORMAT ( 110,5X,F10.3,7X,F10.3)
CONTINUE
GRCUNCWATER FLOW IN FT/DAY
REAC 15,1601 (GwFLOwi i), I»I,NSTEPS>
FORMAT (15F5.2)
WRITS (6,162)
FORMAT (20X, 'GRCUNDWATER FLOW IN FT/OAY')
DO 163 I*1,NSTEPS
WRITE (6,164) I.GWFLOU(I)
FORMAT ( • TIME STEP • , I 5,.'RECHARGE • ,F5. 2)
CONTINUE
REAC (5,195) OELZ.DELT
FORMAT 12F6.2)
WRITE (6,161) DELi.DELT
FORMAT (' SPACING IN CM ISSF8.3 ,5X,.'TIME STEP IN MIN IS',F3.3>
REAC PRESSURE HEADS RATHER THAN MOISTURE CONTENTS
READ (5,190) (THENOTU),J'1,NNODES)
FORMAT (10F6.2)
DTIMEDELT
REAC (5,130) (PSINOT(J),J«1,NNODES)
FORMAT (1PF6.2)
REAC (5,170) (ZTABLEl I ) , I»l .NSTEPS)
FORMAT (1CF6.2)
FOR SANDY SOILS IN STUDY AREA VALUES BELOW ARE ARBITRARY
SINCE TERRAIN IS UNDULATING. ALSQi NO DRAINS EXIST
IN THESE PARTICULAR SANDY SOILS.
READ (5,180) ZDRAIN,ZBOTOM,NOD
FORMAT (2F6.2.I5)
LMGLENOO
READ (5,205) ( HORIZt I ) , I»l, 20)
READ (5,205) ( VERT(I) ,.I »1,6)
FORMAT (20A4)
REAC (5,210) (TIT (I), 1*1,20)
FORMAT (20A4)
REAC IN WITHCRAWAL RATE FOR CROPS IN SANDY SOILS (HEREiCITRUS)
IN I.NCHES/DAY. FROM HASHEMI AND GERBER (1967).
REAC (5,7500) ETCROP.WILTM, DROOT
FORMAT (F5.3t2F5.2>
KEAi:(5,7600) KSCIL
FORMAT (12)
SOIL CLASSES AS FOLLOWS — l.Srt SAND..2. MUCK, 3.NE SAND
CONVERT RAW DATA ON RAINFALLrEVAPORATION,ANO GROUNOWATER FLOW
INTO PROPER AND CONSISTENT UNITS.
CO 200 I=il,NSTEPS
CONVERT INCHES/HOUR TO CM/MIN, BELOW
RAIN(I) = RAIN( I)»0. 0423333
RRAIM I) = RAIN< I )
EVAP(I) = EVAP( I )»0. 0423333
CONVERT FT/DAY TO CM/MIN,.BELOW
GWFLOW(I) * GWFLOW( I)*0. 02116666
CONTINUE
CONVERT CROP WITHDRAWAL FROM IN/DAY TO CM/MIN
ETCROP=ETCROP»0.00176
CONVERT WILTING POINT MOISTURE FROM PERCENT TO FRACTION
WILTM»WILTM/100
• FOR THE 1ST TIME STEP,. COMPUTE VALUES OF VARIABLES
WHICH DEPEND ON PSI
N « 1
NOOSAT(N) = 1
CO 300 J^l.NNODES
PSKN.J) • PSINOT(J)
X » PSINOT(J)
CSLOPE(N.J) = CSL(X)
COND(NiJ) * CK(X)
CONTINUE
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
36
37
38
39
40
41
42
43
44
45
46
47
50
51
52
53
54
55
56
57
58
*Q
3 y
60
61
62
63
64
65
66
67
68
69
70
71
72
73
76
77
78
79
80
263

L = NNODES
CALL VFLUW(N,L,CELZ,ZBOTOM)
LOCATE WATER TABLE
CALL WTABLE(N,L,DELZ»ZBOTOM)
WRITE (6,3000)
WRITE (6,3500) N
CO 2100 J=1,L
WRITE (6,2500) J,PSI(N,J1
2100 CONTINUE
CQ6300J=1»L
X=PSI(N, J)
THETA(N, J) = HUMJCU)
6300 CONTINUE
WRITE OUT MOISTURE VALUES
KRITE(6,5000)
V»RITE(6,5500) N
D06200J=l,L
WRITE(6,6500) J , THETAIN , J >
6200 CONTINUE
PREPARE INFORMATION FOR PLOTTING ROUTINE
INFL
INFL
81
82
NPLOT=L
NSTEP2=L
NCURVE=l
CO 3800 J*1,L
GlKJsj
G(K+1)=TH£TA(N, J)
K=K+2
3800 .CONTINUE
CALL GRAPH 1/2) AND
C SIMILARLY WITH CSLQPt. (RU91N AND STEINHARDT, 1963).
C THUS, FORMULAS FOR PSI ( I ) ,PSI ( 1 1 ) ,P SI U 1 1 I ARE SPECIAL
C CASES OF THOSE USED FOR N 2.
c  „  _  . 
N = 2
NOCSAT(N), * 1
C FIMC VALUES OF PS K I ) , PS I( I I ) , PSI ( 1 1 IJ
CALL PSIMID(NtL)
C FIND COEFFICIENTS A, B, C»D
CALL COEF(N,L,,UELZ,iQOTOM^NOD,.OELT)
C FIND VALUES OF PSI
CALL PSIVAL(N,L,06LT,OELZ)
C FIND VERTICAL VELOCITIES
CALL VFLOW(N,L, CElZ , ZBQTQM)
C LOCATE WATER TACLf:
CALL WTABLE(N,L,DELZ,ZOOTOM)
WRITE (6,3000)
N
2200
WRITE (6,3500)
DO 2200 J=1,L
WR'ITE (6,2500)
CONTINUE
J,PSI(N,J)
X=PSI(N,J)
THETA(N,J)=HUMIC(X)
6<>00 CONTINUE
WRITE OUT MOISTURE VALUES
WRITE16.5000)
WRITQ(6,5500) N
D06600J=1»L
WRITE(6,6500> J,THETA(N»JI
6600 CONTINUE
PREPARE INFORMATION FOR PLOTTING ROUTINE
83
84
85
86
87
88
89
90
91
92
93
95
96
97
98
99
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFL
INFU
INFL
INFL
INFL
INFL
INFL100
INFL101
INFL102
INFL103
INFL104
I NFL 1 05
INFL106
INFL107
INFL108
JNFL109
INFLllO
INFLlll
INFL112
INFL113
INFL1L4
INFL115
INFL1L6
'INFL117
INFL118
INFL119
NPLOT=«L
NST£P2»L
INFL120
INFL121
INFL122
INFL123
INFL124
INFH25
INFL126
IMFC127
INFL128
INFL129
INFL130
INFL131
IMFL132
INFL1J3
264

CO 3900 J = l,L
G(K)=J
G(K+l)=THETA
C FluO VERTICAL VELOCITIES
CALL VFLOWt,N,L,DELZ,ZEOTOM)
C LOCATE WATER TABLE NODE FOR NEXT TIME STEP
CALL WTABLE(N,L,CELZ,ZBOTOM)
WRITE OUT PSI VALUES
WRITE (6,3000)
FOR ALL TIME STEPS
3000 FORMAT f • I1 , 20X , 'TABLE OF PRESSURE HEADS IN CM1)
3500
WRITE (6.3500) N
FORMAT (5X, 'TIME STEP=«, 15)
CO 2 COO J*lfL.
WRITE 16,2500) J,.PS1(N,J)
FOR.VAT ( I5.E8.2)
CONTINUE
HUMICU)
VALUES (PER CENT BY VOLUME)..
2500
2000
C
X=PSI(N, J)
THiiTAIN, J)
4500 CONTINUE'
; WRITE OUT MOISTURE
WRITE (6, 5000),
5003 FORMATC I'iZGX, 'TABLE OF MOISTURE VALUES AS PER CENT BY VOLUME')
WRITE (6,5500) N
5500 FORMAT 15X,'T1ME STEP=»,I5)
U06000J=1,L
WRITE (6,6500) J , TH£TA( N,.J )
FORMAT ( I5.F6.1)
CONTINUE
PREPARE INFORMATION FOR PLOTTING ROUTINE
KKK=l
NPLOTaL
NSTEP2»L
6500
6000
DO 3700 J*l,L
G(K)=.J
G(K+1)»THETA(N,J)
3700 CONTINUE
CALL GRAPH ( NPLOT,.KKK,NSTEP2,NCURVE ,G, TI T)
1000 CONTINUE
CJ.
C WRITE DATA FOR QUALITY PORTION
' WRITE lNUNIF3)DELi»OELT,NSTEPSfNNODES,DRQQT
DO 4700 N=1,NSTEPS
WRITE (NUNIF3)(COND(N,.J)..J*l»NNQDeS)
4700 CONTINUE
CO 4800 N*1,NSTEPS
WRITE (NUNIF3)(VELZ(N,J)rJ"ltNMODES)
INFL134
INFH35
INFL136
INFL137
INFL138
INFLL39
INFL140
INFL141
INFL142
INFLL43
INFL144
INFLU5
INFL146
INFL147
INFL148
INFL149
INFL150
INFL151
INFL152
1NFL153
INFL154
INFL155
INFL156
INFH57
INFL158
INFL159
INFL161
INFH62
INFH63
INFL164
INFH65
TMFI 166
INFH67
INFL168
INFH69
INFL170
INFL171
1NFL172
INFL173
INFL174
INFL175
INFL176
INFL177
INFL178
INFL179
INFU.UO
INFL181
INFL182
INFL183
1NFL184
INFL185
INFL186
INFL187
INFL188
INFL189
INFL190
INFL191
INFL192
265

4800 CONTINUE
CO 4900 N*1,NSTEPS
WRITE (NUNIF3H THETAIN,J),J = 1,NNODES)
4900 CONTINUE
WRITE: (NUNIF3)(QMOLD(N),N=1,NSTEPS)
END FILE NUNIF3
REWIND NUNIF3
RETURN
ENC
SUBROUTINE COEF(NrL»DELZ,.ZBOTOM,NOD,OELT>
C COMPUTES VALUES OF COEFFICIENTS A,B,C,D
COMMON/10 /PSINOT(50),RAIN(20),EVAP(20),GWFLOW(20),ZTA8LE(201 ,
1 PSI<50,20),A(50,20>,B(5C,20).Cl50,.20)fDl50,20)tE(50,.20) »F150,20)
2 CSLOPE(50,20>,THETA(50,20),THENOT(50),VELZ(50,20),PSIH50,20),
3 PSI2(50,.20) ,PSI3(50,20),. NODSATI20),
4 CCNC(50,20),PSIBOT(50),HPONO(50),QSlNK(50,20),QNET12t50)i
5 QPLANT(50,20)
CALL SINK(N,L,NCD,.DPLT,DELZ J
C DEFINE A.B.C.D AT J=l (BOTTOM)
PSI12=(PSI(N,J) +PSHN,J + 11 )/2.
CONC12CMPSI12)
Zl=ZBOTOM*100./3.29
Z2=/.l +OELZ
A(N,J)*1.
Q(N,J)*1.
C(IN, J)=0.
C CROU.NCWATER FLOW IS POSITIVE DOWNWARD
CINiJ)=Z2Zl(GWFLOW(N)*DELZ/COND12)
1  ((QSINK(N,J)»(DELZ/2.)»OELZ)/COND12)
LL = Ll
CO 300 J=.2,LL
PI = PSIKN.J)
P2 = PSI2JN.J)
P3 = PSI3IN.J)
CK1 = CK(P1)
CK2 = CK(P2)
A.(N.J) " CK1/I2»CELZ)
C(N.J) = CK2/(2«DELZ)
Cl = A(N,J)
B2 = ClN^J)
B3 = CSL(P3)»DELZ/DELT
B(N,J) = B1+B2VB3
C VALUES OF QSINK ARE INPUTS TO 0 BELOW
C(N,J) * AIN»J)*PSI(Nlf J*l) * (B3BlB2)»PSUNl,J)
1 * C(N,J)»PSI(N1,.J1) * CK1  CK2  (0 SINK (N, J ) *DELZ)
3(?0 CONTINUE
C DEFINE A.B.C.D AT GROUND SURFACE
A(NiL)=0.
B(N,L)=L.
X=PSI2(N,L)
R=RAIN(N)EVAP(N)
ZLl=ZBOTOM»100./3.28+(L2)*DELZ
ZL=ZL1+DELZ
C NET RAINFALL IS POSITIVE DOWNWARD
C(N,L) = (R»DELZ/CK(X) )ZL*:ZLl  (QSINK(N,Jl»CDELZ/2. )»OELZ/CK(X) )
RETURN
ENO
SUBROUTINE EXTRA IN.U
C USING EXTRAPOLATION FORMULAS. PREDICTS VALUES OF
C PSI FOR THIS TIME STEP.
COMMON/ID /PSINOT(50),RAIN(20),EVAP(20].GKFLOWC20)tZTABLE120)t
1 PSI(50,20),A(50i20>,G(50,20),C(50,20),0(50,20),E<50,20)»F150,20)
2 CSLOPE(50,20>f THETAl 50, 20) .THENOTl 50), VELZ( 50,20) .PSI 1 (50*20) ».
3 PSI2(50,20),PSI3(50,20)r NODSAT(20>,
4 CQNC(59,2Q)tPSIBQT(50),HPONO(5Q}fQSINK(50t20)iQNET12(50}»
266
INFL196
INFL197
INFL198
.INFL200
INFL201
INFL202
INFL203
INFU204
INFL206
INFL207
INFL208
INFL209
INFL210
INFL211
INFL212
INFL2L3
INFL214
INFL215
INFL216
INFL217
INFL218
INFL219
INFL220
IKFL221
INFL222
INFL223
INFL224
INFL225
INFL226
INFL227
INFL228
INFL229
INFL230
INFL231
INFL232
INFL233
INFL234
INFL235
INFL236
INFL237
INFL238
INFL239
INFL240
INFL241
INFL242
INFL243
INFL245
INFL246
EXTR 1
EXTR 2
EXTR 3
.EXTR 5
EXTR 6
EXTR 7
EXTR 8

C
C
C
C
VALUES FROM
FUNCTIONS
PREVIOUS
500
C
C
C
C
C
C
500
5 CPLANTJ50.20)
C0500J=1,L
CELPSI = PSilNl,J)  PSKN2.J)
PSKN,J> = PSHNl.J) * DELPSI
INSTEAD OF EXTRAPOLATING AS ABOVE, USE
TIME STEP TO EVALUATE THE VARIOUS SOIL
PSI(N,J)=iPSI(Nl,J>
CONTINUE
RETURN
END
SUBROUTINE VFLOWfN,L,DELZtZBOTOM)
COPKCN/IO /PSINCT(50),RAIN(20),EVAP(20) »GWFLQW«20»,ZTABtEC20)»
1 PSI(50»20J.A(50»20»,B(50;20UC(50,20),t)t50,20) ,E150,20)>F159,20)
2 CSLQPE150,20),THETA£50,20),THENOT(50),VEL2(50,2Q),PSU(50,20),
3 PSI2(50,20 ),PSI3t50,20)f. NQOSAT(20)»
^ CCJNC(50,20) t.PSIBOTI50),HPONO(50),QSINK(50,20) , ONE T12 150) ,
5 CPLANT(50»20)
APPLICATION OF CARCY'S LAW TO FIND FL£W VELOCITIES.
COMPUTE FLOW VELOCITY BETWEEN LTH AND (LIJTH NODES
THIS FLOW DEFINED AS INFILTRATION
PSITCP= tPS.I(N,LJ+PSI(N.Ll))/2«
X=PSITOP
CONTCP=CK(X)
CCNC(N,LJ=CONTOP
GRTCP=(PSI(N,L)+OELZPSI(N,L1))/OELZ
VELZ(NfU)aCONTOP*GRTOP
RINFtN)=VELZtN,L)
UTILIZATION OF CONDUCTIVITY AT EACH INTERIOR NODE AND TENSION
AT TWO ADJOINING NODES.
COWNWARD VELOCITY DEFINED AS POSITIVE
C0500J=1,L
X=PSItN,J)
COND(N,J)=CK(X)
CONTINUE
NN » L  1
CO 700 J = 2>NN
700
800
900
950
C
C
ZA6V=,(ZBCTOM/3.23J*100.+(Mll»OELZ
PHI
ZBEL=. IZBOTUM/3.28)*100.*(MM1)»DELZ
PHlD£L=ZBEL*PSI(N,Jl)
GRAC=(PHIAa'/PHI3EL )/(2.*DELZ)
VeLZ(N,JI ^ CONCtN, J)»GftAO
CONTINUE
COMPUTE NET FLOW BETWEEN 1ST AND 2ND NODES
J=l
PSU2=(PSI(N,J)*PSHNtJ + in/2.
CGN012=CK(PSI12)
CONO(N,J)»COND12
PHIABV=( ZBOTOM/3.28)«100.«:DELZ + PSI(N,J«1)
PHIBEL=(ZBOTOM/3.28)»100.tPSI(N,JJ
GRACl2=CPHIABVPHI8eLl/OELZ
CNEVl2(fJ)= COND12»GRAD12
VELZ(N,J)*QNET12tN)
WRITE <6,800) N
FORMAT (• l«, IOX ,' DOWNWARD VELOCITIES IN CM/MIN AT TIME STEP",I5)
CO 950 J=5ltL
WRITE 16^900) J,VELZ(N,.J)
FORMAT i I5.E8.2)
CONTINUE
RETURN
END ,
SUBROUTINE WTA8LE1 N,L»DELZ» ZBOTOM)
LOCATES WATER TABLE BASED ON LOCAL FLOW AND
HEAD CONDITIONS
COMMON/ID /PSINCT(50»,RAIN( 20 ), EVAP ( 20) ,GW=LOW<20> , 2TABLE 120) ,
EXTR 9
EXTR 12
EXTR 13
EXTR 14
EXTR 15
EXTR 16
EXTR 17
VFLO 1
,VFLO 3
VFLO 4
VFLO, 5
VFLO 6
VFLO 7
VFLO 8
VFLQ 9
VFLO 10
VFLO 12
VFLO 13
VFLO 14
VFLO 15
VFLO 16
VFLO 17
VFLO 18
VFLO 19
VFLO 20
VFLO 21
VFLO 22
VFLO 23
VFLO 24
VFLO 25
VFLO 26
VFLO 27
VFLO 28
VFLO 29
VFLO 30
VFLO 31
VFLO 32
VFLO 33
VFLO 34
VFLQ 35
VFLO 36
TABE 1
TA8E 2
TABE 3
267

C
C
TABE 2
TABE 21
TABE 22
PSIV 1
PSIV 2
PSIV 3
1 PSI(50,20), A(50,20),B( 50,.2C ) ,.C( 50, 20) ,0( 50,20) ,E (50,20) ,F 150,20) .TABE 5
2 CSLUPE(50,2G),THETA(50,20),THENOT(50),VELZ(50,20>,PSI1<50»20)t TABE 6
3 PSI2<50,20),PS13(50,20), NODSAT(20), TABE 7
4 CCNDt50,20),PSIBOT<50),HPOND(50),QSINK<5C,20),ONEU2(50), TABE 8
5 CPLANT(50,20)  TABE 9
COMMCN/WTAB/JTAB(50),HTAB(50)
C ASSUME THAT NOOi J»l CHOSEN LOW ENOUGH SO THAT IT WILL TABE 10
C ALWAYS BE SATURATED. FREEZE'S EQ.<12) HOLDS, TABE 11
C SCAN PSI VALUES FOR SATURATION FROM BOTTOM NODE UPWARDS. TABE 12
IF (N.EQ.l) HPOND(N)0,
NfJSAT«l TABE 13
DUbOOJ«2,L TABE 14
IF(PSI(N,J).GE.O.) NNSATiNNSAT*! TABE 15
IF(P$l(iJrJ).LT.C.) GO TO 900 TABE^ 16
000 CONTINUE TABE 17
900 NOCSAT(N) » NNSAT
JTABINJ'NODSATtN)
ZTAQLEIN) =ZBOTOM+,
1 PSI<50,20),A150,20),B(50,20),C(50,20),D(50,20),E(50,.20),F(53,20)
2 CSLOPE(50,20),TH6TM50,20),THENOT(50),VELZ(50,20),PSI1(50,20) ,
3 PSI2(50,20),PSI3(50,20)h NODSAT(20),
4 CONC(50,20>,PSIBOT(50),HPOND(50),QSINK(50,20),QNET12(50)*
5 QPLANT(50,20)
COKMON/POND/CPOND(50)
C SET MAXIMUM PONCING DEPTH, IN CENTIMETERS
PCNMAX5.
C TWOPAIN CASES ARISE 1) UPPER NODE UNSATURATED OR R MSAT'D)
C 2) UPPER.NODE SATURATED (WITH A FINITE HEAD)
C CAS5 2 IS STARTED WHEN R SATURATED K,. AND CONTINUES AS LONG AS
C CONTINUITY RELATIONSHIPS INDICATE EXISTENCE OF A FINITE PONDING
C DEPTH AT THE GROUND SURFACE.
IFtfc.EQ.l) HPONC(N)*0.
IF (N.EQ.2) HPGND(N)«0.
C THIS SUBROUTINE GETS CALLED ON 2ND TIME STEP
IF {N.EQ.2).GO TO 700
HPONO(N)HPONClNi) * ((RAIN(N1)EVAP(N1>)VEL2(Nl»L))«OELT
IF (HPONC(N).LT.O.) HPOND(N)«0.
C PONDING DEPTH MUST BE NONNEGATIVE
700 IF(HPGND(N)*GT.3.) GO TO 300
x«o.
CQNSAT'CK(X)
RRAININ)EVAP(N)
C CHECK 'PREDICTED' VALUE OF PSI AT TOP NODE t IF NOT SATURATED
C SCLVE CASE 1.
IF (PSKN.L) .LT.O.) GO TO 850
IF (R.GT.CONSAT) HPOND(N)»(RCONSAT)»DELT
IF(R.GT.CONSAT) GO TO 300
,PSIV 5
PSIV 6
PSIV 7
PSIV 8
PSIV 9
C
C
C
SOLVE CASE 1,BELOW.
850 INOEX'l
NDUPMY»999
MJNSATNCUMMY
CALL MATRIX(N,L, INDEX,NUNSAT,DELZ,)
HPCMD(N),0.
GO TO 150
SOLUTION OF CASE 2, 8ELOW
SET CONDUCTIVITY AT GROUNQ EQUAL TO SATURATION VALUE.
300 IF

PSIGND=HPDNO(N)
PS1GNC=AMAX1(0.,PSIGND)
PSHNtDsPSIGNO
INCEX«2
SEARCH FOR FIRST UNSATURATED NODE FROM TOP
DO 400 1*2, L
K»Lm
IF (PSIU.K).LT.O.) GO TO 500
400 CONTINUE
500 NUNSAT«K
IF INUNSAT.cC.l) NUNSAT»999
CALL MATRIX(N,L,INDEXiNUNSAT,DELZ)
150 DPCNDEI50,201 »F {50,20 1
2 CSLOPEtSO,20)»TM=TA(50i20>,THENOT[50)*VELZ<50,20),PSIl<50»20),
3 PSI2«50,20) »PSI3(5Q,20),. NODSAT(20)t
A CCNC(50»20),PSIBOT(50lHPONO<50>»QSINKt50,20),QNETl2(50)i
5 CPLANT(50,20)
IF (N.EQ.2.) GO TO 510
C04lOJ»ltU
IF U.EQ.l.) GO TO 610
IF (J.EQ.U) GO TO 610
PSIlCN»J)«(PSI(NltJ)«iPSI(NUJimPSI(N»J+mPSI»RAINt20>,EVAP«20)tGWFLOWC20)»*TABLEUOH
1 P$H5Q,20),A<50,20),8( 50f20),C(50i20),D<50,20),E(50t20),Fl50f20)
2 CSLOPE(50,20)tTHETA(50i20),THENOTt50),VEU(50t20),P9IU50«20)f
3 PSI2(50^0)tPSI3(50f20)» NOOSAT(ZO),
4 COND(50,20)tPSlBOT(50)iHPONOt50).QSlNK(50i20itQNET12(SO)«
5 QPLANT(50»20)
COMMON/CROP/ ETCROP. W ILTM
COMMCN/RCOT/CROOT
COMMON/MOLe/CMOLO( 50 I
COMMGN/SOILVKSOU
DEPTH OF SATURATED SECTION A80VE DRAIN
PREVIOUS TIME STEP
C COMPUTE 1ST THE
C FROM RESULTS OF
PMIO
PHID
PHID
PMIO 10
PMID 11
PMID 12
PMIO 13
PMID 14
PMID 15
PflO 16
PMIO 17
PMID 18
PNIO 19
PWIO 20
PMIO 21
PMIO 22
PMIO 23
PMID 24
PMIO 25
PMIO 26
PMID 27
PMID 28
PMID 29
PMIO 30
PMID 31
PMIO 32
PMIO 33
PMID 34
PMIO 35
PMID 36
PMID 37
PMID 38
PMID 39
PMID 40
*$INK 3
SINK 4
SINK 5
SINK 6
SINK 7
SINK 8
SINK 9
269

NOG IS NODE NO AT WHICH CRAIN LOCATED
C
C
C
C
C
C
C
C
IF(PSItNltJ3.GE.O.) KK=KKH
JHPSKNlrJJ.LT.O.) GO TO 600
500 CONTINUE
600 HCRAIN=KK*UELZ/100.
COMPUTE FLOW TO GRAINS (SEE HARRISON)f IN METERS/HOUR.
X=PSI(NlfNQD)
:K(X)»60./100.
C
C
750
700
850
S = 20. 73.28
XLEN=il.
AREA=.XLEN»S
FLOW IN CFS/METER 3ELOW
CCFSPM=QCRAIN«S«XLEN«O.C098
CONVERT TO CFS/FT BELOW
CCFSPF=QCFSPI*/3.28
C
C
SINK 10
SINK 11
SINK 12
SINK 13
SINK 14
SINK 15
SINK 16
SINK 17
SINK IS
SINK 19
SINK 22
SINK 23
SINK 24
CONVERT FROM M/JR TO CM/MIN
GDRAIN=QDRAIN»100./60.
CCNVEKT TO ( 1/MIN)
CDRAJN=QDRAIN/DELZ
NO DRAINS EXIST IN SANDY SOILS OF STUDY AREA
IF (KSQIL.NE.2) QDRAIN=0.
COMPUTE PLANT REQUIREMENTS WITH DEPTH.
TEMPORARILY ASSUMED TO BE ZERO.
COMPUTE PLANT WITHDRAWALS FOR ALL NODES
ASSUME UNIFORM WITHDRAWAL AT TOP HALF OF PROFILE.
WITHDRAWAL AT EACH NODE CANNOT EXCEED AVAILABLE MOISTURE*
AVAILABLE MOISTL'RE = MOISTUKE CONTENT AT NODE (USE VALUE FROM
LAST TIME PfcRiqoi MINUS CONTENT AT WILTING POINT
COMPUTE WITHCRAWALS AT EACH NODE
LRCCT=ORGCT*30.5/D£L2
LLCW=LLRCOT
IF +QDRAIN
SINK UNITS ( L/MINI
CONTINUE
RETURN
END
SUBROUTINE MATRIXIN.L, INDEXrNUNSAT, DELZ )
COMMON/ 10 /PSINCT(50),RAIN(20),EVAP(20),GWFLOW(20),ZTABLEl20)i
PSI(50i20)»AJ50,20),B( 50, 20 ) ,C ( 50, 20 ) ,D( 50,20) ,E ^50^20) .F150.20)
CSLOP£(50,20)tTHETA(50i20),THENOT(50),VELZ(50,20),PSIH50.20) ,
PSI2(50,20) ,PSI3(50,20), NODSAT(20)t
CGNG(50,20) ,PSIBOT(50) , HPQNOl 50) ,QSINK ( 5C, 20) fQNET12(50) t'
CPLANT(50,20)
DIMENSION XLl50i20),RHS(50)
USES GAUSSIAN TECHNIQUE TO COMPUTE PSI VALUES RECURSIVELY*
INITIALIZE LEFT HAND SIDE ARRAY (XL)
CQ550I*ltL
SINK 25
SINK 26
SINK 27
SINK 28
SINK 32
SINK 33
SINK 34
SINK 35
SINK 36
SINK 37
MATR I
, HATR
MATR
MATR
MATR
MATR
HATR
MATR
MATR 10
MATR 11
MATR 12
3
4
5
6
7
8
9
270

XL(1,,J) = 0.
520 CONTINUE
550 CONTINUE
C INSERT VALUES INTO RIGHT H&NO SIDE ARRAY tRHS)
D0600J*1,L
I = J
RHS( I)=0(N,JJ
600 CONTINUE
C INSERT VALUES OF A.BiC INTO LEFT HAND SIDE
C ARRAY (XL)
I»l
JNOCEI
XL(Un«B(N,JNOC£)
XL(I,H1)*A(N,JNODE)
M = L1
CQ660I=.2tM
JNCCE=I
XLI UI1)=C(N, JNODE)
XL(UI)»8(N,JNOC£}
XL(UI+1)=AXL( I«J)/XL( JtJ)
'WRITE 16,1900) JtXL(J,J)
FORMAT (I5,F5.2)
IPREVaI1
JNODE=»I
XL

C «RITE<6,2700) INDEX
2700 FORMAT(2Xi» INDEXES 13)
IF (INDEX. EQ.1J PSI (N, I )=RHS( I J
IF (INDEX. EU. 2) PSICN, I )*PSI(N,.I)
C COMPUTE VALUES OF PSI RECURSIVELY
LL=L1
00570I=.lt.LL
K=Lt+l
C
C
570
C
C
600
610
620
C
C
700
710
720
,315
C
C
500
510
535
520
C
C
JNCUE1=JNODE1
IF ( JNODE1.EQ.NUNSAT) GO TO 580
PSI(NtJNODEl)*RHSUl)XL(Kl,K)»PSKN,JNODE)
GO TO 570
PREVIOUS STATEMENTS APPLYING OARCY'S LAW TO 1ST UNSATURATED NODE
UIRECTLY HAVE BEEN DELETED
CONTINUE
RETURN
END
FUNCTION TEMS(Y)
TENSION VERSUS MOISTURE CONTENT RELATIONSHIPS
FORMULAS BELOW FOR SANDY SOILS
IF (Y.LT.9.5) GC TO 600
IF (Y.GT.42.) GO TO 610
TENS — EXPt (131. 456YI/27. 4112)
GO TO 620
TENS aEXP( ( 13.6879Y)/0.9433)
GO TO 620
TENS EXPM42.5178Y1/0.1335)
IF (Y.GE.43.3) TENS=0.
RETURN
END
FUNCTION CK(X)
CONDUCTIVITY VERSUS PSI RELATIONSHIPS
FORMULAS BELOW FOR SANDY SOILS
IF (X.GE.O.) GO TO 315
C=ADS(X>
IF (G.LT.50.) GO TO 700
IF (G.GT.150.) GO TO 710
CK = 0.5391  0.1056»ALOG(Q)
GO TO 720
CK = 0.1264  0.00009242»ALOG(Q)
GO TO 720
CK = 0.01119  O.COC238i»ALOG(Q)
RETURN
CK=0.126
RETURN
END
FUNCTION CSL(X)
C (DTHETA/UPSI) VERSUS PSI RELATIONSHIPS
FORMULAS BELOW FOR SANDY SOILS
IF U.GH.O.) GC TO 535
P=AQS(X)
IF (P.LT.25.) GO TO 500
IF (P.GT.85.) GO TO 510
CSL = (27.4112/PJ/100.
GO TC 520
IF (P.LE.l.) GO TO 535
CSL = (0.1335/P )/100,
GO TC 520
CSL = (0.9433/P )/100.
GO TO 520
CSL=0.
RETURN
END
FUNCTION HUMID(X)
COMPUTES MOISTURE CONTENT FOR A GIVEN PSI.
FORMULAS BELOW FOR SANDY SOIL.
81
MATR 82
MATR 83
MATR 84
MATR 85
MATR 86
MATR 87
MATR 88
MATR 89
MATR 90
MATR 91
MATR 92
MATR 93
MATR 94
MATR102
MATR103
MATR104
STEN
STEN 2
STEN 3
STEN 4
STEN 5
STEN 6
STEN 7
STEN 8
STEN 9
STEN 10
1
2
3
4
5
6
7
8
9
STEN 11
STEN 12
SCK
SCK
SCK
SCK
SCK
SCK
SCK
SCK
SCK
SCK 10
SCK U
SCK 12
SCK 13
SCK 14
SCK 15
SCK 16
SCSL 1
SCSL
SCSL
SCSL
SCSL
SCSL
SCSL
SCSL 9
SCSL 10
SCSL 12
SCSL 14
SCSL 15
SCSL 16
SCSL 17
SHUM 1
SHUM 2
SHUM 3
272

IF (X.GC.O.) GO TO 330 SHUM 4
S*ABS(X) SHUM 5
IF (S.LT.26.) GO TO 815 SHUM 6
IF tS.GT.86.) GO TO 825 SHUM 7
HUMIC=131. 4527. 41»ALOG

AYB=Yl
CONTINUE
IXA=AXA+.5
IYB=AYB+.5
250 IFUXA^LUO.OR.IXA.GT.IOO) GO TO 260
IF< IYA.LT.O.OR.IYA.GT.50) GO TO 260
CALL PPLOT( IXA, IYA,NSYMiNCT)
260 CONTINUE
YA= (NMAYBAYA) )/< AX8AXA)
IYAaAYA+YA+0.5
IFHXA.LE.IX6) GO TO 250
GO TO 400
C SET PARAMETERS FOR Y DIRECTION
C
290 CONTINUE
IFIAYB.GT.AYAJ GO TO 295
AYB=Yl
AYA=Y2
AXB=X1
AXA=X2
295 CONTINUE
IXA=AXA+.5
IYA=AYA*.5
lYBAYBi.S
300 CONTINUE
IFtlXA.LT.O.OR.IXA.GT.ICO) GO TO 310
IFIIYA.LT.O.OR.IVA.GT.50) GO TO 310
CALL PPLOTCIXA,IYA,NSYM,NCT>
310 CONTINUE
IYA=lYA+\
XA=(N*JAXBAXA) )/(AYBAYA)
IXA=XA*AXA+0.5
N=N+1
IF(IYAIYB) 300i320,AOO
320 IXA * IXB
GO TO 300
400 RETURN
END,
SUBROUTINE CURVE (X, Y,NPT»NCV,NPLOT )
DIMENSION X( 10 1,10), Yl 10 1,10),. NPT( 10 »
COCMON/LAB/ TITLEl20)tXLAB(ll),.YLAB(6),HORlZ(20),VERTI6»
XMAX=,X(li.l)
XMIN=XMAX
YMAX=,YMIN
CO 205 L*l»NCV
IF(NPTM.EG.n) GO TO 205
CO 204 N=il»NPTM
IF(X(N,L).LT.XMIN> XMIN=X(NtL)
IF(X

RANGE=RANGE/(10.««N)
L=RANGE+1.001
206 CONTINUE
IF(L.EG.2) GO TO 209
IFIL.GT.4) GO TC 207
L=4
207 IF

COMMON/LAS/ TITLE ( 20 ) ,XLAB( 11) , YLAB ( 6) , HORIZ ( 20) , VERT 16 J
DATA SVM / 4H**»», 4H+> + +:, 4H'1", 4HXXXX, 4H....f 4H2222i
4H , 4HIIII, 4H  /
IFIK99) 209,220i23Q
200 AI51IY,IX+1)»SYM(K)
RETURN
220 CCNTINUE
I«0
WRITE (6,103) TITLE, NOT
CO 225 11*1,6
1 = 1 + 1
WRITE (6, 101) YL£B( II ),( A< UJ),J=1,101)
IF( II.E0.6) GO TO 228
CO 224 JJ*i,9
II + l
IFd.NE.28) GO TO 221
WRITii(6,i06) VERT(5),VERT(6)r(A(I,J),J»i,101)
GO TO 224
221 CONTINUE
IFd.NE.24) GO TO 222
WRITE (6, 106) VEIU( 1),VERT<2),(A(UJ),J«1,101)
GO TC 224
222 IFl I.NE.26) GO TO 223
WRITE I 6, 106) VERT ( 3 ),VERT(4),{A(UJ)»J» 1,101)
GO TO 224
223 CCNTINUE
WRITE(6,100)IA( UJ),J«1,101I
224 CONTINUE
225 CONTINUE
228 CONTINUE
WR1TE(6,102) XLAD
hR!Tc(6rl05) .HORIZ
100 FORMAH18X» 101A1)
101 FORWAT(F17,3,1X,101A1)
102 FOKMAT(F20.1,10F10.1)
103 FORMAT ( 1H1 , 16X, 20A4* 10X , 'ELEMENT NUMBER1, 14)
105 FORMAT(/40X,10A4)
106 FORMAT(3X,2A4,7X,101A1)
230 CO 250 Ml, 50
CO 240 J*l,101
240 AII, j)»svym
A(I,l)=SYMt8)
250 CONTINUE
CO 260 Jil,101
260 A(51,J)SYM(9J
DO 270 Ml, 101, id
270 A(5UI)SYM(8)
CO 290 Mil, 41, 10
A(1,1)»SYM(9)
AU,101MSYM(9)
290 CONTINUE
RETURN
CND
FUNCTION ACOS(X)
ARC  (1.0  X) / (1.0 + X)
TARG » SQRT(ARG)
ACQS  2.0 * ATAN(TARG)
RETURN
END
SUBROUTINE SOILCU
C PROGRAM ''SOILQUAL1 SLQU 1
c  SLCU 2
C SLCU 3
C PERFORMS WATER QUALITY COMPUTATIONS AT THE SLQU 4
C UNSATURATED SOIL ZONE. PROCESSES CONSIDERED ARE SLQU 5
C ADVECTION, DISPERSION, ABSORPTION AND RELEASE SLCU 6
C CF SUBSTANCES BY THE SOIL. FINITE DIFFERENCE SLQU 7
276

C EQUATIONS UTILIZE PREVIOUSLY COMPUTED VELOCITIES SLCU 8
C AT EVERY NODE IN SOIL PROFILE. AT EACH TIME SLCU 9
C STEPt PROGRAM COMPUTES CONCENTRATION OF SLQU 10
C SUBSTANG6IS) AT EACH NOCE. TRIOIAGONAL MATRIX SLCU 11
C SCHEME IS USED. INTERFACES WITH '1NFILT* PROGRAM. SLCU 12
C • SLGU 13
C PARTIAL LIST OF VARIABLES— SLCU 14
C NNQOESNC. OF NOCES IN SOIL PROFILE SLCU 15
C NSTEPS'WO. OF TIME STEPS IN PERIOD OF STUDY SLCU 16
C CELT^TIME INTERVAL SLQU 17
C DELZ*VERTICAL INTERVAL SPACING SLCU 18
C TSTARTSTARTING TIME OF PERIOD OF ANALYSIS SLCU 19
C TEKD=ENDING TIME OF PERIOD OF ANALYSIS SLCU 20
C FERT(N1«F£RTILIZER INPUT OF N AT NTH TIME STEP SLCU 21
C FERTPIMJFE/UILIZER INPUT OF P AT NTH TIME STEP SLCU 22
C RAINtNJaRAI.NFALL RATE AT N7H TIME STEP SLCU 23
C EVAP«G^OUNCWATER FLOW AT BOTTOM NODE AT NTH TIME STEP SLQU 25
C SQUNAVgRAGf: (TEMPORAL AND SPATIAL) SOIL RELEASE OF NITRATEN
C IN (LB/FTMONTHACRE)
C SCUH»AVERAG£ SOIL RELEASE QF TOTAL INORGANIC PHQSPHOftOUS
C IN (LB/FTMONTHACRE)
C SI\\AV£RAGe SOIL UPTAKE OF NITRATEN
C IN (LB/FTMONTHACRE)
C SlNPAVERAGE SOIL UPTAKE t)F TOTAL INORGANIC PHOSPHOROUS
C IN (LB/FrMONTHACRE)
C SQURMN.DRATE OF N RELEASE BY SOIL AT NTH TIME STEPiLTH NODE SLCU 26
C IN LB/MONTH . $LCU 27
C SCURP(N,L)»RATE CF P RELEASE BY SOIL AT NTH TIME STEPILTH NODE SLCU 28
C IN .LO/MONTH SLQU 29
C SINKN(N,L)»RATE OF N UPTAKE BY SOIL AT NTH TIME STEP.LTH NODE SLCU 30
C IN L3/MONTI: SLCU 31
C SlNKPtN,L)RATE QF P UPTAKE BY SOIL AT NTH TIME STEP.LTH NODE SLCU 32
C IN LB/MONTH SLCU 33
C VELZ^CONC£NTRATION OF MTH SUBSTANCE AT NTH TIME STEP* SLQU 38
C LTH NODE (IN MG/L) SLQU 39
C CKCIN.DDISPERSION COEFFICIENT FOR N»TH TIME ST6Pt LTH NODE SLQU 40
C IN (SO CM/MIN) SLCU 41
C CROPtN)*CROP TYPE GROWN AT NTH TIME STEP SLQU 42
C THETA(NiLIMOISTURE CONTENT AT NTH TIME STEP,LTH NODE SLCU 43
C VISKIN*KINEMATIC VISCOSITY OF WATER SLCU 44
c ALPHASCOEFFICIENT USED TO COMPUTE DISPERSION COEFFICIENT SLCU 45
C (SEE REDOELLtP.16) SLQU 46
C BETA3«COEFFICIENT ALSO USED IN DISPERSION COEFFICIENT FORMULA; SLCU 47
C AS ABOVE SLQU ^8
C CCNC(N,L)»HYDRALLIC CONDUCTIVITY (LENGTH/TIME) FOR SLCU 49
C NTH TIME STEP»LTH NODE SLQU 50
C CCNVFACTOR USED TO CONVERT CONO {CM/MIN) TO
C PERMEABILITY UNITS I SO CM) * 1.91E7 SLQU 52
C CPSRIN.DPERMEABILITY IN UNITS OF UENGTH»»2) SLCU 53
C HEREtSQ CM SLQU 54
C ACIN,.L>»COeF FOR CONCENTRATION AT UPPER NODEt SLCU 55
C NTH TIME STEP AT LTH NODE SLCU 56
C BClNrDCOEF FOR CONCENTRATION AT MIDDLE NODE» SLQU 57
C NTH TIME STEP AT L*TH NODE SLQU 58
C CC(N,.L)»CO£F FOR CONCENTRATION AT LOWER NODE, SLQU 59
C NTH TIME STEP AT LTH NODE SLQU 60
C CCIN.UTERM FOR KNOWN VALUESf SLQU 61
C NTH TIME STEP AT LTH NODE SLQU 62
C BFLUX(N)«FLUX AT BOTTOM NODE AT NTH TIME STEP SLCU 63
C IN LB/AC/YEAR
C CFACTl=CQNV FACTOR FROM LB/MOACRE TO MG/MINACRE10.35
C CFACT2CCNV FACTOR FflOM (OELZ(CM)*1ACRE) TO LITERS4.05«I10»»4) SLQU 66
C CFLCONV FACTOR FROM (CM/tUN)M MG/L) TO ILB/AC/MIN) SLCU 67
277

c
c
c
c
c
c
c
c
c
c
c
c
c
*8.90»(10*»2)
CYEAR=CONV FACTOR FROM (CM/MIN)•(MG/L) TO (LB/ACREYR) * 4.70E03
KMUNTH=MCNTH IN WHICH SIMULATION TAKES PLACE
K=INCEX DENOTING MONTH
CINITIL)=1NITIAL CONCENTRATION AT LTH NODE
FSAT=FRACTION OF SATURATION (PERCENT/100.)
CRCl]T*UEPTH OF ROOT ZONE (FT)
UNIT=UPTAKE RATE OF N BY PLANTS (L8/ACYR)
UPHOS=UPTAKE RATE OF P BY PLANTS (LB/ACYR)
CMCLC(M,N)=CONC. OF MTh POLLUTANT AT MOLE DRAIN AT STEP N
CTA3LE=CONCENTRATION OF SUBSTANCE AT VlATER TABLE NODE
.•_.._._.._••.•..*_»—*—«**p «• — •>—•*•• — •••••'•«*•**»_••*•*• MMV«»M
SET UP COMMON STORAGE AREAS
CCMKCN/IC /NNODES,NSTEPS,DELTrOELZ,TSTART,TEND,FERTN(50),
1 FE«TP150),RAIH150),EVAP(50),GWFLOW(50.) , SOURN ( 50 ,20 ) ,SOURP(50j20)
2 ,SINKN(50, 20),SIMKP(50,20),VELZ(50,20) ,SMOIST,PSI ISO', 20) ,
3 C(2.5C,20),CROP(5G),THETA(50,.20),VISKIN,ALPHA3,BETA3,CIJND(50,20)
4 AC(50,20),BC150,20),CC(50,20),DC(50,20),BFLUX(50),CFL
5,CU4V,K,L,M ,CPER(50,20),CKO(50,20),CFACTl,CFACT2^N,CINITt50)
6 ,CYEAR
CO^MCN/SOIL/KSOIL
COMMCN/CMOLE/CMOLDt 2,50)
COMMCN/LCCMOL/LMOLc
COMMCN/CTAB/CTAQLE(2,50)
COfCCN/OSETS/NGV«IF,NUMIFl,.NilNOVl»NUNREl,NUNGWl,.NUNRE2,NUNIF2»
1 NU^CV2,NUNOV3,NUNOV4,NUNIF3,NUNGW2,NUNGW3,NOVRE1,.NOVRE2,NRERE1,
2 NKER£2,NOVRE3
READ (NUNIF3)OELZ,DELT,NSTEPS,NNOOES,OROOT
CELZ=6.7
CELT=10
NSTEPS=15
NNGDES=10
1,NMONTH)
100
105
130
C
c
c
150
160
165
KMIJNTH=7
READ IN CATA
READ (5,100) (FERTN(K),K=
FORMAT (12F5.2)
READ (5, 100 )(F£RTP(K),.K = 1,NMONTH)
READ (5,105) SOUN»SOUP,SINN,SINP
FORMAT (4F5.2)
READ (5,130) (CROP(K),K=1,NMONTH)
FORMAT (12F5.2)
READ (5,140) VISKIN,ALPHA3,BETA3
,5)
IN SO CM
SMOIST,VMAX,.CFL,CONV,.CFACT1,CFACT2
,5,FlQ.4,F10.9,F10.5,F10.4tF10.2)
(CINIT(L),L=1,NNODES)
.2)
ZEROS, ABOVE, EXCEPT TOP NODE
FORMAT (3F10,
•CCNV« BELOW
READ (5,150)
FORMAT (2F10,
READ (5,160)
FORMAT (10F5,
INITIALIZE AT
.CYEAR
EVENTUALLY MIGHT NEED ONE INITIAL CONCENTRATION VECTOR
FOR EACH SUBSTANCE
READ (5,165) UNIT,UPHOS
FORMAT (2F5.2)
WRITE (6»700) DELZ.DELT
700 FORMAT (• SPACING IN CM' ,.F7.2, • TIME STEP IN MIN SF7.2)
C CONVERT SATURATION MOISTURE VALUE FROM PERCENT TO FRACTION
SMQIST»SMOIST/100.
K'KMCNTH
C CONVERT FROM »PER FT1 TO 'PER NODE' BASIS
C COULD MULTIPLY VALUES BELOW BY SEASONAL AND SPATIAL FACTORS
SLCU 68
SLGU 69
SLGU 70
SLCU 71
SLQU 72
SLCU 73
SL
SLCU 75
SLCU 76
,SLCU 77
SLCU 78
SLQU 79
SLCU 80
SLCU 86
SLCU 89
SLCU 92
SLQU 93
SLCU 94
SLCU 95
SLCU123
SLCU124
SLQU125
SLQU126
SLCU127
SLCU128
SLCU129
SLCU130
SLCU131
SLQU132
SLQU135
278

5/30.5
SCURPCKtL)=SOUP»C£LZ«C.5/30.5
SINKMK,L) = SINN»DELZ»0.5/30.5
SINKPIK,L)*SINP«OELZ«0.5/30.5
NN1=NNODES1
CC 850 L=i2iNNl
SCURP(K,L)=SQUP»D£LZ
SJNKN(K,L)=SINN«DELZ
SIKKP(K,L)=SINP»CELZ
850 CONTINUE
730
/30
/30
/30
185
C
C
C
180
170
C
C
C
C
C
C
C
C
C
c
C«*»
C
210
220
C««*
«»
230
C*»
C«»*
375
SCU*K(K,L)=SOUN«CELZ«0.5/30.5
SI,NKMK,L)SINN«CELZ*0.5/30.5
SINKP(K,L)=SINP*CELZ*0.5/30.5
COMPUTE N AND P UPTAKE BY PLANTS AT NODES WITHIN ROOT ZONE
NRCCT*ORCOT*30.5/DELZ
NLCwER=NNODESNRCOT
IF INLOW6R.LT.1) NLOWER*!
CO 105 L^NLOWER,NNOCHS
IF '( MOWER. EQ10) G0 T0 185
SlNKMK,L»*SINKMK,Ll*CUNIT/12.)*DELZ/tDROOT«30.5l
siNKP(K.L)»SINKP(K,L)t(UPHOS/l2.)*DELZ/(nROOT*30.5)
CONTINUE
00 170 N.liNSTEPS
CO 180 Ll.NNODES
COMPUTE SOURCE AND SINK TERMS FOR EVERY TIME STEP WITHIN
THE CONTh IN QUESTION
RATE UNITS BELOW ARE LB/MOACftE
SOURNIN,L)=SOURN(K,L»
SCURP(N,L)=SCURP(K,L)
SINKN(N,L)=SINKN(K,L)
SINKPIN,L)»SINKP(K,L)
CONTINUE
CONTINUE
FOR ALL TIME STEPS, COMPUTE THE VALUES OF
VARIABLES TO BE USED IN FINITE DIFFERENCE FORMULATION
NOTE—PRESENTLY,. VALUES OF ' SOURN •, • SOURP *, ' SINKN' AND'SINKP*
ARE READ IN. INSTEAD, THESE COULD BE COMPUTED USING FUNCTIONAL
RELATIONSHIP (E.G.,SINUSOIDAL) INVOLVING MIN AND MAX VALUES.
COMPUTE DISPERSION COEFFICIENTS BELOW—
VISKIN=0.00001040 (BASED ON ORLANDO,FLA. MEAN TEMP 71.5 F)
IN (SO FT/SEC.)
CONVERT VISKIN FROM

CO 200 N1,NST£PS
CO 250 L^l.NNODES
C NOTES—
C ALPHA3=38 FOR SAND
C BETA3=1.2 IHOOPES AND HARLEMAN,1965)
C BELCH. CONVERT HYDRAULIC CONDUCTIVITIES TO APPROPRIATE
C UNITS OF PERMEABILITY (FROM CM/MIN TO SQ CM)
CPER(N,L)*=CONO(N,L)*CONV
TERMl±VISKIN»ALPHA3
C CONVERT VELOCITY FROM CM/MIN TO CM/SEC
VZ=VELZ(N,L)/60.
VZ=ABS(VZ)
TERf2=VZ MCP6R(N,L)»»0.5)/VISKIN
IF (TERM2.LE.O.) GO TU 235
C 'CKi;1 BELOW SHOULD HE IN (SC CM/SEC).
CKC(N,L)=TERM1MTERM2)*«BETA3
GO TO 245
235 CKC(N.L)=iO.
C CONVtRT DISPERSION COEF FROM SQ CM/SEC TO SQ CM/MIN
245 CKC(N,L)=CKD(N,L)»6Q.
250 CONTINUE
200 CONTINUE
C PREPARE FINITE DIFFERENCE RELATIONSHIPS FOR EACH TIME STEP
C GOVckMNG EQUATION (ORIGINAL FORM) IS —
C CITHETA»C)/DT = C»D(ThcTA)/C(T) * THETA*D(C)/D(T)
C = V*C(C)/C(Z)
C +CKD(V)*D2(C)/D2(Z)
c + SOUR
C  SINK
C NOTE IN THIS MCCEL. VELOCITY WAS PREDEFINED AS (+) IN DIRECTION
C OF DECREASING NODE NUMBER (DOWNWARD). THUS IN EFFECT THE
C () VELOCITY TERM ABOVE 3ECOMES ( +.).
C PRQCECURE CCLLECT TERMS ASSOCIATED WITH UNKNOWN CONCENTRATIONS
C FOR UPPER, MIDDLE AND LOWER NODES
C CU FOR ALL SUBSTANCES
MPCL2
DO 300 M^J,MPOL
CO 450' N=U,NSTEPS
CALL GROUP
C SOLVE TRIOIAGONAL SCHEME* BELOW
CALL MAT2
450 CONTINUE
C WRITE OUT RESULTS. BELOW
CALL OUTRIT
500 CONTINUE
RETURN
END
SUBROUTINE GROUP
C GROUPS TERMS INTO COEFFICIENTS A,B,C AND D IN TRIDIAGONAL
C SCHEME. UPPER AND LOWER NODES ARE SPECIAL CASES
COMMON/10 /NNQDF,S.NSTEPS,DELT,DELZ,TSTART.TEND,FERTN(50) ,
1 Ft*TP(50).!UIN(5CM, EVAPl 501 .GWFLOW ( 50) . SOURNtSO ,20 ) ,.SOURP (50'.20)
2 ,SINKN(50,20),SINKP(50,20),VELZ(50,20),SMOIST,PSH50,20),
3 C(2,50,20)fCROP(50)«THETA(50,20),VISKIN,ALPHAS,BETA3.CONDI50,20)
4 AC(50,20),BC(50,20),CC(5G,20),DC(50,20) ,BFLUX(5,0) ,CFL
5,CONV,K,L,M ,CPER(50,20),CKD(50,20),CFACTlfCFACT2,.N,CINITt50)
6 .CYEAR
DIMENSION SOURt 50, 20 I ,S INK ( 50,20)
SLCU15T
SLCU158
SLQU163
SLCU164
SLCU165
SLGU166
SLGU167
SLCU16S
SLCU169
SLCU170
SLQU171
SLCU172
SICU173
SLCU174
SLOU175
SLCU177
SLCU178
SLQU179
SLCU180
SLCU181
SLCU182
SLCU183
3LCU184
SLCU185
SLQU186
SLCU187
SLCU188
SLCU189
SLCU192
GRUP
GRUP
GRUP
SL
GRUP
GRUP
.GRUP
GRUP
GRUP
1
2
3
5
6
7
8
9
C
c
c
300
150
CIRCUMVENT CALCULATIONS FOR 1ST TIME STEP
INSTEAD, INITIALIZE CONCENTRATION VALUES
IF (N.GT.l) GO TO 150
CO 300 L^l.NNODES
C(M,N.LHUNIT(L)
CONTINUE
GO TO 250
SET SOURCES AND SINK VALUES ACCORDING TU SUBSTANCE
DO 350 LiliNNODES
GRUP 10
GRUP 11
GRUP 12
GRUP 13
GRUP 14
GRUP 15
GRUP 16
GRUP 17
280

SUBSTANCE 1
IF (M.EQ.l)
IF (N.EO.l)
SUBSTANCE 2
IF (M.EQ.2)
IF (M.EQ.2)
350 CONTINUE
IS NITRATE NITROGEN
SOUR(N,L)=SOURN(N,L)
SINK(N,L)=S1NKN(N,L)
IS TOTAL INORGANICPHOSPHOROUS
SOUR(N,L)=SOURP(N,L)
SINK/OELT)
* (THETA(N,L)/DELT) +: IVELZ(N,L)/OELZ)
'CC(NL>=0.
CC(N,.L)= (THETA < N, L ) *C (M,.N1,L)/DEL T>
L +.( (SOUR ( Ni DS INK(N,L) >*CFACT1 /1CFACT2* (DELZ/2. ) »THETA(N,L)) 1
C
C
PERFCftM FOR INTERMECIATE NODES
CU 100 L=i2iLL
AC(NiL)=(VELZ(N,L)/(2.»DELZ) )+.( CKOlN.L ) / IDELZ»»2) )
3C(NiL>=( (THETA(NtL)THETA(Nl,L),)/DELT) + ( THETA(N»L) /DEIT)
I + (2.»CKOIN>L)/(DELZ»»2) )
CC(NtL)(VSLZ(N,L)/«2.*DELZ))*(CKD(NL)/(DELZ#*2))
CC(N,L)= (THETA(N,L)»C(M,.N1,L)/DELT)
1 t ( (SOUR (N,L)SINK ( N, L ).) »CFACT1/ (CFACT2«OELZ»THETAIN, L))
100 CONTINUE
C
C
C
C
c
c
c
PERFORM FOR TOP NODE ,L=NNODES
L=NSCDES
AC(NiL)=0.
BC(N,L>= «ThETA(N,L)T»ETA(Nl,L))/DELT) * (THETA (N,L) /CELT)
I  (VELZtN.D/OELZ)
CC(NrL)»  VELZIN.D/OELZ
CC(N,L)= (THETA(N,L^•C(M^N1,L)/DELT)
1 *{(SOUR(NtL)SINK(N,L))«CFACT1/(CFACT2*(DELZ/2.)»THETAtNiL)»)
nMfSINKP(50,20),VELZ(50,20) , SMOl ST, PSH 50,20) ,
3 C(2,50t20),CROP(501tTHETA<50,20),VI SKIN,ALPHA3,BeTA3,COND(50(20)
4 AC«50t20),BCC50,20),CC(50,20>,DC(50,20).BFLUXl50)iCFL
5,CONV,K,L,M ,CPER(50,20),CKO(50,20),CFACTlhCFACT2,NtCINITt50)
6 ,CY£AR
CO^VCN/LOCMOL/LfOLE
CCMMC'J/CMCLE/CMCLO( 2i 50)
COMKCN/WTAB/JTABI50)HTAB(50)
CO^CN/CTAB/CTA6LE( 2i50)
DIMENSION XL(5Ci20»,RHS(50)
CIRCUMVENT CALCULATIONS FOR 1ST TIME STEP
IF (N.£(,».!> GO TO 530
USES GAUSSIAN TECHNIQUE TO COMPUTE CONCENTRATION VALUES
INITIALIZE LEFT HAND SIDE ARRAY (XL)
C0553I=1,NNODES
C0520J = liNNODES
XL(UJ) = 0.
520 CONTINUE
550 CONTINUE
C
C
INSERT VALUES INTO RIGHT HAND SIDE ARRAY (RHS)
C0600J=1,NNODES
I*J.
GRUP 18
GRUP 19
GRUP 20
GRUP 21
GRUP 22
GRUP 23
GRUP 24
GRUP 25
GRUP 27
GRUP 28
GRUP 29
GRUP 30
GRUP 31
GRUP 32
GRUP 33
GRUP 35
GRUP 37
GRUP 38
GRUP 39
GRUP 40
GRUP 41
GRUP 42
GRUP 43
GRUP 44
GRUP 45
GRUP 47
GRUP 48
GRUP 49
MTUO
MTWO
MTWQ
MTWO
SL
MTWO
MTWO
,MTWO
MTWQ
1
2
3
4
6
7
a
9
MTWO 10
MTWO 11
MTWO 12
MTWO 13
MTWO 14
MTWO 15
MTWO 16
MTWO 17
MTWO 18
HTWO 19
MTWO 20
MTWO 21
MTWO 22
MTWQ 23
MTWO 24
281

600
RHS( I)=DC(N,J)
CONTINUE
C
c
c
660
C
1400
C
C
200C
CISCO
INSERT VALUES OF A,B,C INTO LEFT HAND SIDE
ARRAY (XL)
1=1
JNCCE*!
XL«I,I)=BCM
JNODE=»I
XL»

570
C(P,N,JNODEl)
CONTINUE
'RHS(Kl)  XL(K1,.K)»C(M,N,JNOOE)
C
C
MTWO 90
MTWO 91
MTWQ 92
C
C
C
SET CONC. AT MOLE DRAIN BELOW
LMQLE=i
CMOLD(M,N)=CIM,N,LMOLE)
SET CONC. AT WATER TABLE NOCE
JTAB(N)=1
JJ=JTAB

C •*••••••«*»» •••••••^••••» • ••^ J» «• *»W ^ffc^«**^it»^» <••**• ^«* ^••••i
C CONSIDER OUTPUT TO OVERLAND FLOW MODEL
C 
WRITE (NUNOV4>NTIME,DTIME
WRITE (6..110) NTIME,OTIME
110 FORMAT (I5.F5.1)
C CONVERT RRAIN FROM CM/MIN TO IN/HR
DO 100 N^ltNTIME
RRAIN(N)=RRAIN(N)*60./2.5*
WRITE (NUNOV4)RRAIN(N)
WRITE (6,120) RRAIN(N)
120 FORMAT (10X.E8.2)
100 CONTINUE
C CCNVERT RINF FROM CM/MIN TO CM/HR
CO 200 NMtNTIME
RINF(N)»RINFIN)«60.
WRITS (NUNOVORINFH)
WRITE (6,120) RINF(N)
200 CONTINUE
WRITE (NUNOV,N«liNTIME)
WRITS (6»130) (DPONO(N) ,N*1,NTIME)
130 FORMAT (5X.10E8.2)
END FILE NUNOV4
REWIND NUNOV4
C CONSICER CUTPUT TO RECEIVING MODEL IN MUCK FARMS—
C 
C FLUW AND QUALITY IN MOLE DRAINS
IF (KSOIL.NE.2) GO TO 300
WRITE (NUNRE1) NTIME.OFIME
C WRITE FLCW IN CFS PER UNIT FOOT OF MOLE LENGTH
WRITE (NUNRE1) ( QMOLOJ N ) ,N*i,NTIME )
C WRITE NFLUX IN L8/DAY PER FT OF MOLE LENGTH
M»l
CO 400 N*1,NTIME
C CONVERT MGD»MG/L TO LB/OAY
FLUXN(N)«.{QMQLO(N)/1.5S)*(CMOLD(H,N) )»8.34
WRITE (NUNRE1) FLUXN(N)
400 CONTINUE
C WRITE PFLUX IN LB/DAY PER FT OF MOL6 LENGTH
K»2
DO 500 N»sl,NTIME
FLUXP(N)^(QMOLC(N)/l.SSt*(CMOLD(M,N))*8.3A
WRITE (NUNRE1I FLUXP(N)
500 CONTINUE
C TO COMPUTE FLOWS AND FLUXES FOR EACH CANAL JUNCTION,
C yUST MULTIPLY ABOVE VALUES BY (AREA OF CATCHMENT(SQ FT) /20t)
C SINCE SPACING OF MOLE DRAINS IS APPROXIMATELY 20 FT.
300 CONTINUE
C CONSIDER OUTPUT TO GROUNDWATER MODEL
C 
C WATER TABLE ELEVATION AND CONCENTRATIONS—
C TRANSFER VALUES AT LAST TIME STEP
NN»NTIME
CMNN)*CTABLEUNN)
CP(NN)"CTABLE(2,NN)
C WRITE WATER TABLE HEIGHT AND NITROGEN CONC. ONLY
WRITE (NUNGW3) IT3
WRITE (NUNGW3) CN3
END FILE NUNGM3
REWINC NUNGW3
RETURN
ENU
SUBROUTINE OVLAN1
COMPON/DSETS/riGWlF,NUNIFl,NUNOVl,NUN*El,NlJNGWlt.NUNRE2,NUNIF2V
1 NUNOV2,NUNaV3,NUNOV4,NLiNIF3,NUNGW2,NUNGW3,NOVRELf.NOVRE2,NREREl,
284

c
c
c
c
c
c
c
c
c
c
c
c
c.
c.
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
2 NRERE2.NOVRE3
RETURN
END
SUBROUTINE OVLAN2
OVERLAND FLOW PROGRAM »KWAVE'
THIS MODEL USES KINEMATIC WAVE RELATIONSHIPS
TO COMPUTE FLOWS FROM PLANE CATCHMENTS TO
RECEIVING ELEMENTS. CATCHMENTS 00 NOT FEED INTO EACH OTHER.
A RECEIVING ELEMENT MAY BE CONNECTED TO MORE THAN ONE
CATCHMENT. SCANNING PROCEDURE USED TO OBTAIN
EQUIVALENT STORM DURATION, AVERAGING PROCEDURE
US£0 TO OBTAIN EQUIVALENT 'MEAN (CONSTANT) RAINFALL
ANC INFILTRATION RATES. QUALITY RELATIONSHIPS ARE
INCLUDED. PROGRAM APPLIED TO SANDY SOILS IN STUDY AREA
(APCPKA BASIN).
PARTIAL LIST OF VAR1ABLES
NRbC»NO. OF RECEIVING ELEMENTS (TOTAL)
NDEP.NO. UF RECEIVING ELEMENTS CONSTITUTING «SMALL»
SURFACE DEPRESSIONS (TYPED
NfiCCY NG. OF RECEIVING ELEMENTS LOCATED ON
•LARGE' BOCY OF WATER (TYPE2)
NCATCHUJNO. OF CATCHMENTS FEEDING TO ITH
RECEIVING ELEMENT
NTYPEIIJTYPE OF ITH RECEIVING ELEMENT (1 OR 2 ABOVE)
XU)*XCOOROINATE OF ITH RECEIVING ELEMENT
Y(M«YCGQRDINATE OF ITH RECEIVING ELEMENT
NTCTOTAL NO. OF CATCHMENTS
XLMIK,M)'LENGTH THROUGH KTH INTERNAL IN KTH CATCHMENT
•IN FEET
XLENGJKHLENGTH OF KTH CATCHMENT I IN THE HORIZONTAL)
IN FEET
WIOTH(K)»WIDTH CF KTH CATCHMENT
IN FEET
SLOPE(K)«SLOPE CF KTH CATCHMENT
•IN (FT/FT)
ALPHA(K)«ALPHA FACTOR OF KTH CATCHMENT (SEE EAGLESON)
(OlfENSIONLESS)
EMK)»MANNING 'N' ROUGHNESS COEF FOR KTH CATCHMENT
EM(K)«MFACTOR OF KTH CATCHMENT (SEE EAGLESON)
(HEREi EMIKJ2.)
NLIN'KJKJa NO. OF RECEIVING ELEMENT TO WHICH KTH CATCHMENT
IS LINKED
MINTING. OF INTERVALS INTO WHICH EVERY CATCHMENT IS DIVIDED
NT»NO. OF TIME STEPS IN PERIOD
CTTIME INCREMENT IN MINUTES
XRAIN(L)«RAINFALL RATE AT LTH TIME STEP
IN (INCHES/HOUR)
X1NFIDINFILTRATION RATE AT LTH TIME STEP (FLOW BETWEEN UPPER
TWO NODES IN UNSATURATED MODEL) IN (CM/HOUR)
XPCNCIDPONOING DEPTH AT LTH TIME STEP (FROM UNSATURATEO
MOCEDt IN CM
LSTOPTIME STEP AT WHICH PONDING DEPTH BECOMES ZERO
XFAVERAGE INFILTRATION RATE FOR I L LSTOP
IN (FT/MIN)
NTRNO. OF TIME STEP AT WHICH MAXIMUM PONDING DEPTH OCCURS
TRTIME CORRESPONDING TO ?NTR' ABOVE*ASSUMED TO BE EQUIVALENT
STORM DURATION, IN MIN
XIxAVERAGI? RAINFALL RATE FOR 1 L LSTOP
IN (PT/MIN)
YlENGm^ELEVATION DIFFERENCE BETWEEN UPPER ANDLOWER PORTIONS
OF KTH CATCHMENT,. IN FT
SL(K)»»TRUE« LENGTH OF KTH CATCHMENT
IN FEET
SLM(K,M)*TRUE LENGTH,AS ABOVE.THRU MTH INTERVAL
IN FEET
285
1
2
3
4
5
6
7
8
9
KINW
KINW
KINW
KINW
KINW
KINW
KINW
KINW
KINW
KINW 10
KINW 11
KINW 12
KINW 13
KINW 14
KINW 15
KINW 16
KINW 17
KINW 18
KINW 19
KINW 20
KINW 21
KINW 22
KINW 23
KINW 24
KINW 25
KINW 26
KINW 27
KINW 28
KINW 29
KINW 30
KINW 31
KINW 32
KINW 33
KINW 34
KINW 35
KINW 36
KINW 37
KINW 38
KINW 39
KINW 40
KINW 41
KINW 42
KINW 43
KINW 44
KINW 45
KINW 46
KINW 47
KINW 48
KINW 49
KINW 50
KINW 51
KINW 52
KINW 54
KINW 55
KINW 56
KINW 57
KINW 58
KINW 59
KINW 60
KINW 61
KINW 62
KINW 63
KINW 64

C SINTH(K)SINE QF ANGLE THETA WITH HORIZONTAL FOR KTH CATCHMENT KINW 65
C TC(KrM)=TIME OF CONCENTRATION IN KTH CATCHMENT KINH 66
C AT MTH POINT, IN MIN KINW 67
C VEL(K,L.M)=OUTFLOW VELOCITY IN KTH CATCHMENT AT OUTLET KINW 68
C AT LTH TIME STEP AND AT MTH POINT, IN FT/MIN KINW 69
C YL(K,.L,M)=OEPTH IN KTH CATCHMENT AT LTH TIME STEP,AT MTH POINT KINW 70
C IN INCHES KINW 71
C CUNIT(K,L,M)=UNIT OUTFLOW,KTH CATCHMENT AT LTH TIME STEP, KINW 72
C AT MTH POINT* IN (CUBIC FEET PER MIN/FT) KINW 73
C CTCT,(K,L,.M) = CUTFLOW AT KTH CATCHMENT AT LTH TIME STEP, KINW 74
C AT MTH POINT, IN (CUBIC FT / SEC) KINW 75
C CRCV(I,L)<=INFLOW INTO ITH RECEIVING ELEMENT AT LTH TIME STEP KINW 76
C IN (CUBIC FE6T/SECI KINW 77
C TI(K,M)=T1M£ AT WHICH YL VANISHES IN KTH CATCHMENT.AT MTH POINT KINW 78
C IN MIN KINW 79
C TP(K,M)=TIME AT WHICH YL BEGINS TO DECREASE AFTER ACHIEVING KINW 80
C PEAK DEPTH, AT KTH POINT, IN MIN KINW 81
C VBAR(K,L)=AVERAGE VELOCITY AT KTH CATCHMENT AT LTH TIME STEP KIKW 82
C IN FT/SEC KINW 83
C YBAR(K,L)=AVERAGE DEPTH AT KTH CATCHMENT AT LTH TIME STEP,IN FT KINW 8*
C PMS(K,L)=RATE OF INPUT (SUSPENSION) OF SOLIDS INTO WATER KINW 85
C IN LB/ACRE/SEC KINW 86
C PMD(K,L)=.RATE OF INPUT OF DISSOLVED SUBSTANCES INTO WATER KINW 87
C IN LB/ACRE/SEC KINW 88
C R=HYCRAULIC RADIUS, IN FT KINW 89
C VCRIT=AVERA3E VELOCITY AT WHICH PARTICLES OF DIAMETER *0« KINW 90
C WILL BE TRANSPORTED, IN FT/SEC KINW 91
C CK=SHIELC'S CONSTANTS.C56 KINW 92
C SS=SPECIFIC GRAVITY OF PARTICLE (APPROX.=2.) KINW 93
C C=CIAMETER OF PARTICLE (ASSUMED TO Be D50) ,IN FT KINW 94
c PS«AX=MAXIMUM SCOURING RATE IN LB/ACRE/SEC KINW 95
C CPOL(K,L)*FLUX OF SUSPENDED POLLUTANT (PHOSPHOROUS IN THIS CASE)
C IN OUTLET FLCW (L3/OAY)
C RO=DENSITY OF WATER KINW 98
C PZESC=INITLAL RATE OF UPTAKE OF DISSOLVED SUBSTANCE (NITRATEN
C IN THIS CASE) IN ,T1(200,20UTPI200,1),XLM(200,I),MINT,SLMt200 ,1),
6 YL(200,20,l),VEL(200,.20,l), QTOT(2CO,20U) »
7VBAR(200,20),YBAR(200,.20),PMS(200,20),PKD(200,20),CKrSS,D^PSMAX,
8 CPOL(2()0,20).RO,PZEROtCOECAY,CDIS(200,20)
9 ,X(200),YI200),I,K,L,M
COK^CN/DSETS/NGWlF,NuNIFl,NL.NOVl,NUNREl,NUNGWlrNUNRE2,NUNIF2i
1 NUNCV2,NUNOV3,NUNOV4,NUNIF3».NUNGW2,NUNGW3,NOVREl,NOVRE2iNRERElt
2 NRERE2fNOVRE3
C,
C
C
C
C
READ INPUT DATA BELOW
REAC PHYSICAL DATA AND LINKAGE DATA BETWEEN RECEIVING ELEMENTS
AND CATCHMENTS
286
KINW116
KINW1X7
KINW118
KINWH9
KINW120

C
C
C
C
C
C
100
120
125
110
135
130
140
C
C
C
150
160
350
450
210
170
180
READ (5,100) NREC,NOEP,NBODY,.MINT,NTC
FORMAT (515)
CO 125 I*1,NREC
READ (5,120) II,NTYPEm,NCATCH(I),X(I),Y(I)
FORMAT ( I3,I2,12,4X,2F5.2)
CONTINUE
READ (5,110) (NCATCHII),1=1,NREC)
FORMAT (2512)
DO 130 K*1,NTC
REAC (5,135) KK,NLlNK(K),.XLENG{K)^WIOTH{K)rSLOPE(K)
. ,EN(K)(EM(K)
FORMAT (14, 3X,I3,2F5.0,F5.4,F5.3,F5.2)
SET VALUES OF EMK) FOR TURBULENT FLOW
£M(K)=5./3.
KC(K)sKK
CONTINUE
REAC HYDROLOGIC DATA
READ (5,140) NT.CT
FORMAT ( I5,F5.2)
FOR THE PRESENT, ASSUME OATA BELOW HOLDS FOR EVERY CATCHMENT
DO 160 Lnl,MT
READ (5,150) XRAIN(L),XINF(L),.XPONO(L)
FACTOR*!.3
ADJUST RAINFALL INTENSITY FOR VARIOUS RUNS
XRAIN(LI»XRAIML»*FACTO.R
UNITS FOR 'XRAIMIN/HR ,. FOR 'XINF 'CM/HR , FOR «XPOND'*CM
NOTE— VARIABLES ABOVE COULE BE OBTAINED THROUGH INTERFACING
INSTEAD OF READING IN
FORMAT (3F10.2)
CONTINUE
REAC INUNOV4)NT,DT
CO 350 L*1»NT
PEAO (NUNOV4) XRAIN(L)
CCUTINUE
WRITE (6»2IO) (XRAINtL),L*l,NT)
CO 450 L*1»NT
Rfc'AC (NUNOV4) XINF(L)
CONTINUE
WRITE (6,210) (XINFU),L*1,,NT)
READ (NUNOV4)(XPONO(L).L*1,NT)
WRITE (6,210) {XPONDtL),L«l,NT)
FORMAT (5X.12E8.2)
REWIND NUNOV4
REAC (5,170) CK,SS,0
FORMAT (2F5«l,F7.5l
READ 15,180) PSMAX,PZERO,CDECAY
FORMAT (3F5.1)
CALL AVE
WAVE
RECV
RQUAL
WRIT
CALL
CALL
CALL
CALL
RETURN
C DEBUG SUDCHK,TRACE,INIT,SUBTRACE
END
SU&rtOUTlNE AVE
C PERFORMS AVERAGING COMPUTATIONS INVOLVING
C RAINFALL AND INFILTRATION. PREPARES DATA FOR USE 'IN
C KINEMATIC WAVE CALCULATIONS.
COMMON/ID / NR£C».NDEP,NBODYrNCATCH(70>,KC(200J
1,NTYPE (2 ),NTC,XLENG( 200),.WIDTH! 200)» SLOPE 1200) t
2 iLPHA(200),EN(200),EM(200),NLPJK(200»,NT,DT,XRAIN(200),
3 XINF(2C»,XPONO<20>,ISTOP,XF,NTR,.TR,XI,
4 YLENG(200),SLt200),SINTh(2CO),TC(200t1).
5 QRCV(200,23)tTI(200,20),.TP(200,l),XLM(200,l)«MINT,SLMI200»l)»
6 YL<200,20,lJ,VEH200i20,.l), QTOT(20020,ll f
7VBAR(200r20),YBAR(200,20),PHS(2COf.20),PMO{200,20)
8 CPCL<200*20),RC»PZERO,CDECAY,CDIS(200».20)
287
KINH121
KINW122
KINW123
KINW124
KINW125
KINW126
KINW127
KINM128
KINW129
KINW130
KINW131
KINW132
KINWL33
KINW134
KINW136
KINW137
KINW138
KINM139
KINM140
KINW141
KINW142
KINW1A3
KINW144
KINW147
KINW148
KUW149
KINW150
KINW151
KINH152
KINM153
KINW155
KIKW1S6
AVE 1
AVE 2
AVE 3
AVE t

9 iX(200) ,.Y(200),ItK,L,M
C SCAN TIME HORIZON TO FIND TIME STEP AT WHICH AVE 15
C PONDING CEPTh BECOMES 'PERMANENTLY' ZERO. WORK AVE 16
C PACKWARD IN TIKE. AVE 17
CG 200 L=a,NT AVE 18
NDUNKYsNT+1L AVE 19
IF fKC(200l
ltNTYPE(2)fNTC,XLENG<2CO)t.WICTH<2QO),SLOPE{200),
2 ALPHAl200).EN(200)iEM(2CO)tNL!NK(200),NT,DTiXRAIN<200)»
3 XINF120),XPONC(20),LSTOP,XF,,NTR,TR,XI,
4 YLENG(200)tSL(200),SlNTH{2CO»,TC(2C0.1)r
5 ORCV(200,20),TI(200,20)i.TP{200tl),.XLM(200tl)tMINTfSLMl200ill*
6 YL«200,20?l),VEL(200i2C».l) , QTOT (200,20 ,1) »
7VBAR(200»2C) »YBAR( 2CO,20),PMS( 2CO,20),PMD(200,20) ^C K ,SS ,D t PSMAX ,
ti CPCL(20J,20>,RG,PZEROrCDECAY,CDIS(200t20)
9 fX(200)r.Y<200)t IiKiLtH
CO 200 KnltNTC WAVE 17
YLEhGlK)=XLcNG(K)*SLOPE(K) WAVE 18
RSC*XLENG(K)»»2 4 YLENG(K)»»2 WAVE 19
SL

C TURBULENT FLOW ASSUMED ABOVE
200 CONTINUE
PERFORM THE FOLLOWING CALCULATIONS FOR EVERY INTERVAL
CO 700 M=U,MINT
COMPUTE TIME OF CONCENTRATION BELOW
CO 300 K^ltNTC
SLK(K,M)=SL(K)*(M/FLCAT(MINTJ)
XNUM=SLM(K,M)« UXIXF)»M l.EM(K)) )
CONVERSION BELOW— 60 SEC/MIN
XC£N=ALPHA(K J »60.
EMINV=1.0/EM(K)
TC(K,M)= (XNUM/XOEN) »*EMlNV
300 CONTINUE
WRITE (6,3501 TR
350 FORMAT CIS1 STORM DURATION*' ,F7. 2, 'MIN* )
CO 360 K=a,NTC
WRITE (6,370) K,TC(K,M)
370 FORMAT (' CATCH N0= ' , 15 , • TIME OF CONC=',F7.2,
360 CONTINUE
C FOR EACH CATCHMENT COMPUTE DEPTH AT OUTLET
C
C
C
DEPENDING ON RELATIVE VALUES OF TR AND TC(K).
CASE 1 TR TCIK,M)
CASE 2 TR TC(K,M)
CO 400 K^ltNTC
IF (TR.GE.TC(K,M))
IF (TR.LT.TCtK.M))
SUBROUTINE HIDUR
SUBROUTINE LOOUR
CALL HIOUR
CALL LOOUR
*MIN» )
400 CONTINUE
C CCNPUTE VELOCITIES AND FLOWS USING DEPTHS
CO 500 K=U,NTC
DO 6CO L=a,MT
C 'VEL1 BELOW IN FT/SEC
C (COMPATIBLE KITK MANNING'S N IN ALPHA COMPUTATION*
C PREVENT BASE ZERO IN OPERATION BELOW
If (YL(K>L,M).LE.O.) YL (K ,.L ,M) =0. 1E15
VbL(K,L,M) = ALPHA(K.) » I YL( K ,L,.M ) »• ( EM(K )L. > )
C CLNVERT TO (FT/MN)
VeL(K,L,M)=VEL(K,L,M)»60.
C 'CUNIT' BELOW IN (CU FT PER MIN/FT)
CUNIT =VEL(K,L,M) « YL(K,L,,M)
C 'OTOT' BELOW IN (CU FT/SEC)
C3»0.0167
QTCT(K,L,M)=CUNIT • WIDTH(K)*C3
600 CCNTINUE
500 CONTINUE
700 CONTINUE
RETURN
C DEBUG SUBCHKtTRACE,IN1T.SUBTRACE
END
SUBROUTINE RECV
C COMPUTES TOTAL FLOW INTO EACH RECEIVING ELEMENT
C FROM EACH CATCHMENT,FOR EACH TIME STEP
COMtfCN/ID / NREC,NCEP,NBODY,NCATCH{70),KC(200)
L,NTYP£<2),NTC,XLENG(2CO),WIQTH<2QO),SLOPE<200)t
2 ALPHA(200),EN(200),EM(200),NLINK(200)tNTtOVrXRAINI200) t
3 XINF(20),XPCND{20),LSTCP,XF,NTR,TR,XI,
4 YLSNG(200),SL(200),SINTH(2CO),TC(200,l)l
5 QRCV(200,20),TI( 20(J, 20),.TP( 200,1 ),XLM(2CC,1) ,MI NT, SLM(200il) ,
6 YL(200,20,l),VEL(200,20t.l)t QTOT (200,20,1) t
7VBA^(200^20),YBAR(200,2C>,PMS(200,.20},PMO(200»20)»CK>5S,U»PSHAX»
8 CPCL(200,20),RC,PZERO,CCECAY,COIS(200,20)
9 ,X(200)^Y(200),I,K,L,M
C WORK ONLY WITH LAST INTERVAL OF EACH CATCHMENT
M*MINT
C INITIALIZE ARRAY BELOW TO ZERO
CO 200 L=»1,.NT
CO 300 I*1,NREC
289
WAVE 23
WAVE 24
HAVE 25
WAVE 26
HAVE 27
WAVE 28
WAVE 30
WAVE 30A
WAVE 31
WAVE 32
WAVE 33
WAVE 34
WAVE 35
WAVE 36
WAVE 37
WAVE 38
WAVE 39
WAVE 40
WAVE 41
WAVE 42
WAVE 43
WAVE 44
WAVE 45
WAVE 46
WAVE 47
WAVE 48
WAVE 50
WAVE 51
(HAVE 53
WAVE 54
WAVE 55
WAVE 56
HAVE 57
RECV 1
RECV 2
RECV 3
RECV 14
RECV 15
RECV 16
RECV 17
RECV 18

CRCV(I,L)=0.
300 CONTINUE
200 CONTINUE
CO 400 L=1,NT
CO 500 K»itNTC
INDEX=NLINK(K)
•CJRCV BELOW IN ICU FT/SEC)
CRCVlINDEX,L)=QRCV(INDEXED
500 CONTINUE
400 CCNTINUE
RETURN
END
SUBROUTINE HICUR
UTILIZES KINEMATIC
TR TC(K,y) [HIGH
+ QTOT(KtL»M)
FOR THE CASE
C UTILIZES KINEMATIC WAVE RELATIONSHIPS
C TR TC(K,y) [HIGH CURATION)
CCMKCN/IC / NREC,NOEP,N30DY,NCATCH(70),KC(200)
1,NTYPE(2J,NTCtXLENG( 200 ) ,W ICTH( 200) .SLOPE (200) ,
2 ALPHA(200), EN(200),EM(200),NLINK(200),NT,DT,XRAINt2001 .
3 XI;F(20),XPCND(20),LSTOP,XF,NTR,TR,XI,
4 YLEiiG(200),SL(2CO),SINTH(2GO)tTC(200tl) ,
5 (JRCV(2QOt20),TI<200i20),.TP(200,l),XLM<20Ci 1) .MINT , SLM (200 .1) t
6 YLC200,20,l),VEL(20Q,20,.l). OTOTl200.20,1)t
7VBA3(200,20),YBAR<200,20),.PMS(200,.20),PMD(200,20) ,CK*SS,D«PSMAX,
8 CPCL(200,20),RC,PZERO,COECAY,CDIS(200,20)
9 ,X(200) »Y(200)t I,.KTL,M
C CETrRMINE THK,M),TIME AT WHICH YUK.L.M) VANISHES
THK,M)«TR +

c
c
c
CASE 1 OELOW
COMPUTE DEPTH AT DIFFERENT TIMES
100 CO 150 L^1,.NT
T=L*CT
110 IF IT.GT.TR) GO TO 120
120
130
140
150
GO TC 150
IF (T.GT.TP(KtM)) GO TO 130
YL{KfL,M)»XI*TRXF«T
GO TO 150
IF (T.GT.THK.M) ) GO TO 140
TERM1=(XIXF)«»0.5
TERM2=(XI»( *0.
CONTINUE
RETURN
• (XIXF))*«0.5
C
c
c
CASE 2 BELOW
COMPUTE DEPTH AT DIFFERENT TIMES
200 CO 250 L*l•NT
T=L»CT
210 IF (T.GT.TR) GC TO 220
YL(K,L,M)=(XIXF)*T
CO TO 250
220 IFU.GT.TPCK.M) ) GO TO 230
YL(K,LiM)teXI*TRXF«T
GO TC 250
230 YL(K,L,M)=0.
250 CONTINUE
RETURN
C
C
C
CASE 3 BELOW
COMPUTE DEPTH AT DIFFERENT TIMES
300 CO 350 L^liNT
T=L»CT
310 IF (T.GT.TR) GC TO 320
YL(KrLtM)*tXIXF)»T
GO TO 350
320 IF (T.GT.THKiM) ) GO TO 330
YL(K*L,M)=XI»TRXF*T
GO TC 350
330 YLCKtL,M)=0.
350 CONTINUE
RETURN
END
SUBROUTINE RCUAL
C PERFORMS WATER QUALITY COMPUTATIONS AT THE VARIOUS
C SUBCATCHMENTS. RELATES THESE TO INFLOWS INTO
C RECEIVING ELEMENTS.
COMMON/IC / NRECiNOEPfNBODYtNCATCHt70)»KC(200)
1 ,NTYPE(2 ) ,NTCtXLENG( 2CO ) ,WlDTH( 200)» SLOPE (200) *
2 ALPHA! 200), EN( 200), EMI 2CO),NLlNM2CO)t NT tDTtXRAlNt 2001 t
3 XINF(20),XPCNO(20),LSTOP,XF,NTR,TR,XI,
4 YLENGI200)tSL(2CO),SlNTH(2CO>,TC(200tl),
5 QRCVt200.20),TI<200,20)iTPC200,l),XLM(200,U,MINT,SLMI200»lU
6 YL(2CO,20, i),VEL(200t20.l), QTOT<200,20t1 ) »
:7VBAR(200t20>.YBAR(2COt20)fPMS(200f20},PMD(200»20),CK,SS,D^PSMAX,
a CPULI200.20J,RO,PZ6RO,CDECAY,CDIS(200i20)
9 ,X«2GOUY<200) , ItK.LtM
COKMGN/NUM3/QPHCS(2CO,20) ,QNIT(200,20)
C
C
C
AT EACH TIME STEP,
CATCHMENT BASED ON
DO 100 L*liNT
COMPUTE MEAN VELOCITY FOR EACH
VELOCITIES AT ITS INTERVAL POINTS.
291
LOOV 23
LQOV 24
LOOV 25
LCOV 26
LODV 27
LCOV 28
LODV 29
LCDV 30
LCDV 31
LODV 32
LOOV 33
LOOV 34
LODV 35
LCOV 36
LOOV 37
LCDV 38
LCDV 39
LODV 40
LODV 41
LCCV 42
LCDV 43
LODV 44
LOOV 45
LCCV 46
LCDV 47
LODV 48
LODV 49
LCDV 50
LCDV 51
LCOV 52
LODV 53
LOOV 54
LCCV 55
LODV 56
LODV 57
LCOV 50
LCDV 59
LODV 60
LCDV 61
LGDV 62
LODV 63
LCDV 64
LODV 65
LODV 66
LCDV 67
LCDV 68
LODV 69
LODV 70
LCDV 71
RQUA 1
RQUA 2
RCUA 3
RCUA A
RCUA 15
AQUA 16
RQUA 17
RCUA 18

DO 200 K=il,NTC
M
C
c
C
c
c
c
c
c
c
c
c
c
c
c
c
SUM2»0.
CO 300 M*il, MINT
SUM2*SUM2*YL(K,L,M)
300 CONTINUE
C40.0167
•V8AR' BELOW IN FT/M1N
VBAR(K,L)»(SUM1/MINT)
•YOAR1 BELOW IN FT
YBArUK,L)«SUM2/MINT
20»0 CCNTINUE
100 CONTINUE
COMPUTE VALUES CF PMS(K,L) AND PMD

570
TADJ^T UXLENG(K)/2.)/VBAR/43560.
CSFLUX UNITS IN (LB/SEC)
CCIS(K,L>=DSFLUX«86«00.
IF OUTFLOW RATE IS ZERO » SET COIS»0.
IF (QTOT2^)fTI(20020)^TP(200tl)»XLM(200l}fMINT*SLM(200Ut*
YL(200,20»l>fVEU(20Q,2C,l), QTOT(200(20il)t
7V&AR(200»20>tYBAR(200,2G),PMS(200i20)iPMO(200t20}»CKfSStO»PSMAXi
8 CPCL(200t20)tROtPZERO,CDECAY».CDIS(200,20)
9 ,X(200>F.Y(20P),IiK,L,M
COM^ON/NUM3/CPHCS(200i20) »GNITI200*20)
M»MINT
WRITE (6,300)
FORMAT (' STEP
WRITE (6t350)
FORMAT (•
CO 100 K^l.NTC
CO 200 L*ltNT
WRITE (6^250) LtKiQTOT(Kf.L*.M)*CPOL (K,L) »CDIS(K*L)
FORMAT (2I5,3F10.2)
CCNTINUE
CONTINUE
RETURN
END
SUBROUTINE OVLAN3
COMMON/DSETS/NGKIF,NUNIFltNUNOVUNUNREl*NUNGWl»NUNRE2*NUNIF2t
1 NUMOV2tNUNOV3iNUNOV4»NUNIF3*NUNGW2,NUNGW3,NOVREl,NOVRE2*NRERElt
2 NRERE2.NOVRE3
RETURN
END.
293
RQUA 68
RCUA 69
RCUA 70
RCUA 71
RGUA 72
RQUA 74
RQUA 75
RQUA 75B
ftCUA 77
RCUA 78
RCUA 79
RQUA BO
RGUA 81
RCUA 82
RGUA 83
WRIT 1
WRIT 2
5
6
CAT FLOW
P LOAD
N03N LOAO»J
WRIT 13
WRIT 14
(CFS) (LO/DAY) (LB/OAY) ')
WRIT 18
WRIT 19
WRIT 20
WRIT 21
WRIT 22

SUBROUTINE OVLAN'.
CUKlrfON/DSETS/NGV

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c..,
PHIPRV = PREVIOUS VALUE OF PHI, IN FT.
PHINEW = UPDATEC VALUE OF Phi, IN FT.
PHICCR = CORRECTED VALUE OF PHIt IN FT. , AFTER RELAXATION
DELPHI =• DIFFERENCE BETWEEN PRESENT AND PREVIOUS VALUE OF PHI
TOL = ABSOLUTE VALUE OF DELPHI, ABOVE (TQLERENCE)
FRAC=FRACTION OF NODES FOR WHICH AN ERROR LARGER
THAN TOLERANCE VALUE IS ALLOWED
NPER = NO. OF NCCES CORRESPONDING TO FRACTION
PHISYM = SPECIAL VALUE OF PHI DETERMINED FROM
ABOVE
SYMMETRY FORMULA
PHILO = LOWER BCUNO VALUE USED IN EQUIPQTENTI AL CALCULATIONS
PHIHI = UPPER UCUND VALUE OF PHI USED IN EOUl
POTENTIAL CALCULATIONS
REMUREM2 » REMAINDER TERMS USED TO DETERMINE LOCATION OF EQUI POTENTIALS
XORIG * VALUE OF XCCOKDINATE AT ORIGIN
OHPHI = HORIZONTAL GRADIENT OF PHI
DVPHI = VERTICAL GRADIENT OF PHI
VELH = HORIZONTAL VELOCITY { FT/DAY )
VELV = VERTICAL VELOCITY < FT/DAY )
VELR = MAGNITUDE OF VELOCITY RESULTANT
VH,VV,VR  ABSOLUTE VALUES OF VELH, VELV, VELR
SINA=SIN£ OF ANGLE OF RESULTANT VELOCITY
THEFA*ANGLE OF RESULTANT WITH RESPECT TO POSI
ABOVE
TIVE
HORIZONTAL AXIS, COUNTERCLOCKWISE, IN DEGREES
ANGLE ~ MATRIX USED TO STORE VALUES OF THETA ABOVE
IDIV = NO. OF DIVISIONS ALONG HORIZONTAL
JDI7 = NC. OF DIVISIONS ALONG VERTICAL
XMI.N » MINIMUM VALUE OF ABSCISSA ,IN FT.
XtfAX * MAXIMUM VALUE OF ABSCISSA , IN FT.
YMIN = MINIMUM VALUE OF ORDINATE , IN FT.
YMAX = MAXIMUM VALUE OF ORDINATE , IN FT.
L = fATRIX USED TO STORE VALUES OF KCLASS » ABOVE
U = HORIZONTAL VELOCITY (POSITIVE TO RIGHT)
VY ~ VERTICAL VELOCITY (POSITIVE DOWNWARD)
DCL '. DISPERSION COEF ALONG DIRECTION OF FLOW
CCT * DISPERSION CUEF TRANSVERSE TO DIRECTION
OF FLOW
NUMBERING SYSTEM FOR DISPERSION NODE CLASS MATRIX KCLASS
0 = EXTERIOR
I * LEFT GOUNCARY
2 » LEFT TOP
3 =» LEFT BOTTOM
4 = TOP
5 = DCTTOM
6 = RIGHT TOP
7 = RIGHT BOTTOM
8 = RIGHT BOUNDARY
9 = INTERIOR
10 = KNOWN CONCENTRATION
PHI 7
COK»«CN/ZEROTH/C,CSTAR,CSS., I I, JJ , OT, OX, DY,DCL,DCT,L,
IA1,A2,A3,.A

c
c
c
c
c
c
c
c
c
c
c
c
c
READ (NUNGW2). HT2
HT2=C.
REAC (NUNGW3). HT3
HT3=0.
READ IN NITROGEN CONCENTRATIONS
REAC (MUNGWl) CM
CN1=0.
PEAC (NUNGW2) CN2
CN2=l.
REAC (NUMGW3) CN3
CN3=0.
REWIND NUNGWI
REWIND NUNGW2
REWINC NUNG/J3
REAC (5,100) NROW,NCOL
WRITE (6, ICO) NROW,NCOL
100 FORMAT (213)
•NCLASS' CODE IS 0 0 FOR
2 FOR AOOVE WATER TABLE,
HEAD.
INTERIOR
3 FOR NO
PHI 20
PHI 21
PHI 22
NODES,I FOR WATER TABLE* PHI 23
FLOW BOUNDARY, 4 FOR KNOWN PHI 24
PHI 25
CO 110 J=a,NROW PH4 26
C PEAC IN FROM BOTTOM TO TOP
REAC (5,120) (NCLASS(I,J> ,I»l,NCOL) PHI 27
120 FORMAT (2011) PHI 29
110 CCNTINUE PHI 30
C WRITE (6,127)
127 FORMAT (• FLOW CLASS MATRIX")
CO 125 J=1,NROW
JJ=NROWJ+1
C WRITE (6*126) ( NCLASSU , JJ ), 1 = 1 tNCOL)
126 FOkMAT (2011)
125 CONTINUE
CO 133 I=1,NCOL PHI 31
F6AD (5,140) (CROUNCU UWTABLEU),BOTOM(1,I),BOTOM(2,I)) PHI 32
C WRITE (6*141) (GROUND( I ),.WTABLEU),BQTOK.EQ.l) PHItI,J)=PHI{I»J)+HT1
181 CCNTINUE
296

CO 182 1=^13,14
CG 182 J5.UNROW
IF (NCLASS.C I,J).EQ.1» PH II I, J >*PHI( I, J )+HT2
182 CCNTINUE
CO 183 I=a5,NCQL
CO 183 J=1,NROW
IF INCLASS(IiJ).EQ.l) PHI(I,J>PHI(I,J)+HT3
183 CONTINUE
REAC (5,200) OELX,D£LZ,S,ZZERO
C WRITE (6,200) DELX,OELZ,S,ZZERO
200 FORMAT (4F10.2)
REAC (5,210) ITMAX.TOLER,OMEGA
C WRITE (6,211) ITKAX,TOLER',OMEGA
210 FORMAT (15,2F5.4)
211 FORMAT (IX,I5.2F5.2)
C SPECIFY ELEVATION OF LOWER BOUNDARY,BELOW.
READ (5,220) ELECT
C WRITE (6,220) ELBOT
220 FORMAT (F6.1)
PEAC (5,230) NREG
C WRITE (6,230)NREG
230 FORMAT ( ID
C READ J COORDINATES FOR EACH REGION,STARTING WITH
C LOrtER LEFTHAND CORNER,THEN CLOCKWISE.
CU 240 1=1,NREG
REAC (5, 250) UK I),K1( I ),L2(I),K2(I ) ,L3 < I ) ,K3( I ) ,LM I) ,K4l!) )
C WRIT£(6,250)(LUI),K1(I),L2(I),K2( I),L3(I),K3(I),L4(I),K4LI) )
250 FQRPAT (812)
240 CONTINUE
CO 270 1^1,NREG
REAC (5,260) ( PHI ( LU I ) ,K1( I) ), PHI (L2 ( I) ,K2( I ) ).
I PhI(L3(I),K3(I>), PHI(L4(I),K4(I)))
C WRITE (6,260) ( PHI(LllI),K1(I)), PHI(L2(1),K2t PHI(L4(I),K4(I)» )
260 FORMAT (4F6.2)
270 CONTINUE '
REAC15.2BO) ( JT ABLE( I), I*1,.NCOL )
CO 290 JROW=1,NROW
REAC (5,285) (KCLASS(I,JROW),I=1,NCOL)
C READ IN FROM TOP TO BOTTOM IN THIS CASE
285 FORfAT (2012)
290 CONTINUE
280 FORMAT (1912)
C.
C
COMPUTE VALUES THAT HILL REMAIN PERMANENT
ALPHA = (CELX/DELZ)»»2
COMPUTE HEAD FOR ALL NODES IN EACH RECTANGULAR REGION
CO 300 M=1,NREG
PHIREG*tPHI(Ull*>t.KllM) )«PHI(L2(M)»K2(K) )+PHI(L3(M) ,K3lM)»
+PHI1L4
PHI(I,J)=PHIREG
CONTINUE
CCNTINUE
CONTINUE
GO TO 320
GO TU 320
GO TO 320
WRITE OUT INITIAL VALUES OF PHI
WRITE (6,325)
325 FORMAT (30X,•INITIAL PHI MATRIX")
CO 330 J»ltNROW
K=NRCH*1J
PHI 61
PHI 62
PHC 63
Phi 64
PHI 65
PHI 66
PHI 67
PHI 68
PHI 69
PHI 70
PHI 71
PHI 72
PHI 73
PHI 74
PHI 75
PHI 76
PHI 77
PHI 78
PHI 79
PHI 80
PHI 81
PHI 82
PHI 83
PHI 84
PHI 85
PHI 86
PHI 87
PHI 88
PHI 89
PHI 90
PHI 91
PHI 92
PHI 93
PHI 94
PHI 95
PHI 106
PHI107
PHU08
PH1109
PHI110
PHI111
PHI112
PHI113
PH1114
PHI11S
PHU16
PHI117
PHI118
PHI119
PHI120
PHI121
PHI122
PHI123
PHI124
PHI125
PHI126
PH1127
PHI128
PHI129
PHI130
297

c.
c
WRITE (6,340) (PHUIiK) , I=1,NCOL)
340 FORMAT (20F5.1)
330 CONTINUE
"*"sTART*ITERATIVE PROCESS TO FIND PHI VALUES
NNCDES=NCOL«NRQW
CO 400 M=5l,ITMAX
c
c
GO TO 700
GO TO 700
C
C
C
C
C
C
C
540
550
700
500
450
C
C
C
c
c
c,
c
850
800
400
CO 450
CO 500 I*1,NCOL
SCAN TO CETERMINE NOCE CLASS
PHIPRV=PHI
IF (\CLASS(I,J).EQ.4)
TYPc BELOW IS IMTERIOR
IF (NCLASS

980 FORMAT (5X..20I5)
CO 990 J=1,NROW
995
990
WRITE (6,995» ( Ji ( PhK UK ) , 1 =1 ,NCOL ) )
FORMAT ( I5,20F5.11
CONTINUE
CALL EQUIP( I,J,C£LX,CELZ,ELBOT,NROW,NCOL)
C INSERT 'CLEANUP1 ROUTINES PRIOR TO PLOTTING
C INSERT PLOTTING ROUTINES
CALL VEL( I,J,DELX,OELZ, NROW.NCOL)
CALL PRECIS ( I, J, DELX.DELZ.NROW.NCOL.ZZERQ)
RETURN
END
SUBROUTINE BOUNC < I, J, ALPHA)
C THIS SUBROUTINE SELECTS THE PROPER NUMERICAL EXPRESSION
C TO BE USED FUR ANY 'BOUNDARY* NODE.
COKMCN/FINIT/NCL*SS<20,20),PERKH120,20) , PERKV(20t 20) ,PHI(20,20> ,
1 JTABLE(20),KCLASS( 20,20)
C POINTS ON BOUNDARY ARE THOSE SATISFYING I=1,I=NCOL
C AND/CR J=l. NOTE J=l HAS 2 CORNER POINTS.
IF. (J.NE.l) GO TO 1000
C FORMULAS BELOW FGR Jt
IF (l.EQ.l) GO TO 1200
IF (I.EQ.NCOL) GO TO 1300
GO TO 1400
12CO Xi\UMER= (2»PERKH I, J )»PHIl 1 + 1, J ) *2*ALPHA»PERKV< I , J) »PHI ( H Jf 1) )
XCcKCK= 12*PERKH I,J)+2«ALPHA»PERKV(I,J))
GO TU 1500
1300 XNUMER=(2»PERKH( I1,J)»PHI< 11 , J ) t,2»ALPHA«PERKV( I » J ) *PHI (IiJ+l>.>
XDENCM= (2«PERKH IliJ) +! 2*ALPHA*PERK VU ,J ) )
GO TO 1500
1400 XNUMER= (PERKHt Il»J)«PHUlfJ)+PERKHUiJ)»PHI

4000
C
C
C
5500
5000
C
C
PHIUiJ)* (PHI(IltJ)4PHI( !4ltJ})/2.
GO TO 5000
IF IHHKI*liJ).EQ.O.J GO TO 5500
PERFORM FOR LEFT NODE PKI»0.
PHISYH=IPhIt IfJ*l) + PHU I1,J1) ) /2.
FORMULA BELOW GOOD FOR HOMOGENEOUS MEDIUM ONLY.
XNUMER = (8»PHISYM44»PHI( Hl.J) 4. 3» (PHI ( I *J«1) +
1 PHIt I,J1I I 1
XDcNCJM»18.
PHK I,J)=iXNUMER/XDENOM
GO TC 5000
PERFORM FOR RIGHT NODE PHI'C.
PHISYM  IPHl(I,J4.1)4PHlU4;UJl)}/2.
XNUMER» (8*PHISYM+4*PHI( I1,J) * 3* ( PHI ( I , J*i) +
1 PHI )
XCENOM»lfl.
PHId, JlXNUMER/XDENOM
RETURN
ENCr
SUBROUTINE EQUIP ( I , J,DELX, DELZ ,ELBOT,NROU,NCOL)
THIS SUBROUTINE COMPUTES ELEVATIONS (MSL) OF
POINTS OF KNOWN HEAD, 6Y INTERPOLATION, FOR EACH COLUMN
TABL 13
TABL 14
TABL 15
TABL 16
TABL 17
TABL 18
TABL 19
TABL 20
TABL 21
TABL 22
TABL 23
TABL 24
TABL 25
TABL 26
TABL'27
TAOL 28
TABL 29
TABL 30
TABL 31
EGUI 1
EGUI 2
EQUI 3
COMMCN/FINIT/NCLASS( 20, 20 ) , PERKH( 20 » 20 > »PERKV( 20 ,20 ) ,.PHI{20»20) , ECUI 4
C
C
C
105
C
C
400
300
200
100
450
C
1 JTABLEI20),KCLASS(20,20)
SPECIFY LCWER AND UPPER BOUNDS FOR HEAD, BELOW
PHILO»50.
NL*PHILO
PHIHI125.
NHaPHlHI
SCAN ALL COLUMNS FOR EACH HEAD VALUE
HEKE, INCREMENT IS 5 FT.
WRITE (6,105)
FORMAT (•!',//, 5X, 'LOCATION OF EQUIPOTENTIAL POINTS IN
1 ' )
CO 100 LlNL,NH,5
DO 200 Ul.NCOL
JTAB'JTABLEI IJ1
CO 300 J^liJTAB
REMUPHldi J)L
IF (REM1.GE.O.) MM»l
IF (REM1.LT.O.) PM«0
REM2=PHId»J + l)L
IF (REM2.GE.O.) NN»1
IF (REM2.LT.O.) NN*0
KC ( 1 U v lU U A. hi h4
\&Vi>it^'\^t*
'IF (NSUM.NE.l) GO TO 300
WH£N HEAD VALUE FALLS BETWEEN NODES, INTERPOLATE
TO FIND ELEVATION
RI»ABS(REMi)
R2»ABS(REM2)
ZINTn (Rl«DELZ)/(Rl+R2)
ZL» ELBOT 4 U1)«OELZ 4 ZINT
IF (NCLASS (I.JJ.EQ.2) GO TO 300
IF(NCLASS(I, J41).£Q.2) GO TO 300
WRlTt*(6, 400 ) L,I»J,ZL
FORMAT (2Xt »PHI»S I3» 'COL"1 , 13, 'ROW*' , 1 3t>*ELEV«* ,F7.2)
CONTINUE
CONTINUE
CONTINUE
WRITE (6,450)
FORMAT ( '1', 20X, 'INTERPOLATIONS IN THE HORIZONTAL')
SCAN ALL ROWS FOR EACH HEAD VALUE
XQRIG'O.
ITMAXNCOLl
CO 500 L*?NL,NHf 5
CO 600 JMtNROW
CO 700 IM,ITMAX
REM.la*PHI ( It J )L
ECUI 5
EGUI 6
EQUI 7
EGUI 8
EGUI 9
EGUI 10
EQUI 11
EQUI 12
EGUI 13
GRID SYSTEMEGUI 14
EQUI 15
EQUI 16
EGUI 17
EGUI 18
EQUI 19
EQUI 20
EGUI 21
EGUI 22
EGUI 23
EQUI 24
ECUI 25
EQUI 26
EGUI 27
EQUI 28
EGUI 29
EGUI 30
EGUI 31
EQUI 32
EGUI 33
ECUI 34
EGUI 35
EQUI 36
EQUI 37
EQUI 38
EQUI 39
EQUI 40
EGUI 41
EGUI 42
EGUI 43
EQUI 44
EGUI 45
EGUI 46
EQUI 47
EQUI 48
EQUI 49
300

c
c
INTERPOLATE
800
700
600
500
C
C
50
IF (REM1.GE.O. MM»l
IF (REMl.LT.O. MM.0
REM2PHI I I+:i,J L
IF (REM2.CE.O. NN«1
IF IREM2.LT.O. NN0
NSUMHM*NN
IF (NSUM.NE.l) GO TO 700
WHEN HEAD VALUE FALLS BETWEEN NODES,
TO FIND HORIZONTAL POSITION.
Rl»ABSIREMl>
R2ADSIREM2)
XINTMR1«OELX)/(R1*R2)
XLXCRIG+:i I1)*CELX+XINT
IF INCLASSiI,J).EQ.2) GO TO 700
IF (NCLASS11+1,J).EQ.2) GO TO 700
WRITE (6,800) LiI.J.XL
FORMAT (2X,.«PHI»,.I3,ICOLI,I3» 'ROW*1,1 3r«X*» »F8.1)
CONTINUE
CONTINUE
CONTINUE
RETURN
END
SUBROUTINE VEL ( I,J,CELX^DELZ.NROW.NCOL)
THIS SUBROUTINE COMPUTES THE VELOCITY VECTOR AT
EACH NODE IN THE GR 1C SYSTEM (EXCEPT AT THE BOUNDARY NODESk)
COKJ"!ON/PINIT/NCLAS$(20t20}fPERKH(20,20),PERKV(20i20)»PHI<20,20)
1 JTAELE(20),KCLASS120,20)
COMMGN/VELOZ/VELH(20,20)».VELV<2C!,20UVELR(20»20),ANGLE120»'20)
COPMCN/ZERaTh/C,C$TARfCSS,Il,JJ,DT,.DX»DYf.DCL*DCT,L,
lAl,A2,A3,A4,BltB2».83,KFX,.CFX,KLX>CLX,KFY,CFY1KLY,CLY»
2U»VY,ZPiGPiGPC»
3XMIN,XMAX,IDIV,YMIN,YMAX,JOIV
DIMENSION CI20,20 ) ,CSTAR(20,20),CSS(20,20»«U(20,20),VYl20*20)
DIMENSION L(20,20)
INITIALIZE VELOCITIES TO ZERO
CO 50 I»1»NCCL
DO '50 Jl.NROW
VELH(I,JJ»0.
VELV(I,.J)»0.
CONTINUE
DO 100 tl.NCOL
CO 200 Jil.NROW
'STAR* OUT ALL NONINTERIOR NODES.
IF (NCLASS(I,J).EQ*2) VELVlI,J)«9$99999.
VELHU,J>«9999999.
VELRl UJ) "9999999.
ANGL£(I,J)9959999.
GO TO 200
(NCLASSIItJ).E0.2)
(NCLASSt I»J).ECJ.2}
(NCLASSIIiJ).EQ.2>
[NCLASSd, JJ.EQ.2.)
IF
IF
IF
IFI
SET HORIZONTAL AND VERTICAL VELOCITIES AT ZERO AT NOFLOW
IF (NCLASSI I,J).EQ.3) VEIVU,J)»0.
EGUI
EOUI
EQUI
EQUI
EGUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
EQUI
ECUI
EQUI
EQUI
EGUI
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
IF
IF
IF
IF
(NCLASSII,J).EQ.3)
(NCLASSIIiJ).EQ.3)
(NCLASSII,J).EU.3)
(NCLASSII,J).EQ,3)
C
C
150
VELH(I,J)*0.
VELR(I,J}0.
ANGLEU, J)»0»
GO TO 200'
JRQWNROW+1J
IF (KCLASS(I,JRCW).NE.9> GO TO 150
USc 2POINT FORMULA FOR NUMERICAL DIFFERENTIATION
FIND PHI GRADIENT IN HORIZ..VERTICAL
CHPhI — PHI(IfJm/f2.*DELZ)
COMPUTE VELOCITIES IN FT/DAY
VELHII,J) • PERKHII,J)*0.134*DHPHI
VELV(I,J) * PERKV(I,J)*0.134*DVPHI
GO TO 180
DETERMINE HORIZONTAL GRADIENTS
IF (KCLASSII,JRCW).EQ.2) DH0.
IF iKCLASStI,JROW).EQ.A) DH0.
IF (KCLASSII»JROW).EQ.6) DH0.
80UNDARYVEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
VEL
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
1
2
3
4
5
6
8
9
10
11
12
12A
12B
12C
120
12E
12F
13
14
15
16
17
18
19
20
20 A
20B
20C
200
20E
20F
21
22
23
24
25
26
27
28
29
30
301

c
c
IF
IF
IF
IF
IF
FOR
(KCLASSt I,JROW).ECU1)
(KCLASSt I» JROWJ..EQ.3)
(KCLASSII,JRCW).EQ.7)
(KCLASS1I,JRCW).EQ.8)
(KCLASS(ItJROW).6a.5)
DH=(PHI< I+:i,J)  PHI (UJn/DELX
DH*(PHI
KCLASSI It JRCWJ.EQ.4)
KCLASSt I,JROW).EQ.6)
FOR POLLUTANT INPUT POIN
CALCULATIONS IN 'GWQUAL1
300
400
450
700
500
600
200
100
800
850
910
VH=VcLH(I,JJ
VH= ABS(VH)
VV = VELV( UJ)
VV= AES(VV)
VR = VELR(I»J)
VR = ABS(VR)
IF (VR.LT.0.1E12) VR=Q.1E12
SINA = VV/VR
A = ARSIN(SINA)»180./3,1416
START SEPARATING INTO QUADRANTS
NDUNMY = MNSUM1.
IF(NCUMMY) 300,400t500
THETA = 180. + A
3RD CiUADRANT
GO TO 600
IF(VELV< ItJI.LT.O.) GO TO 450
CISTINGUISH BETWEEN 2ND AND 4TH QUADRANTt ABOVE
THETA = 180.  A
GO TC 700
THETA i 360.  A
GO TC 600
THETA = A
1ST GUADRANT
ANGLEtItJ) = THETA
CONTINUE
CONTINUE
WRITE OUT VELOCITY AND ANGLE
WRITE (6,800)
FORMAT (•!•,//,30X,'VELOCITY FT/DAY  AND ANGLE  DEC  MATR1XM
WRITE (6,850) tIfI=l,NCOL)
FORMAT (2016)
CO 900 J=1,NROW
K=NRCW+1J
WRITE (6..910) tVELR(l,K),l = l,NCOL)
FORMAT (20F6.U
WRITE (6,920) (ANGLE!I,K),I=1,NCOL)
VEL 34
VEL 35
VEL 36
VEL 37
VEL 38
VGL 39
VEL 40
VEL 41
VEL 42
VEL 43
VEL 44
VEL 45
VEL 46
VEL 47
VEL 48
VEL 49
VEL 50
VEL 51
VEL 52
VEL 53
VEL 54
VEL 56
VEL 57
VEL 58
VEL 59
VEL 60
VEL 61
VEL 62
VEL 63
VEL 64
VEL 65
VEL 66
VEL 67
VEL 68
VEL 69
VEL 70
VEL 71
VEL 72
VEL 73
VEL 74
VEL 75
VEL 76
VEL 77
VEL 78
VEL 79
VEL 80
VEL 81
VEL 82
VEL 83
VEL 84
302

920 FORMAT (20F6.1,/M
900 CONTINUE
RETURN
ENU
SUBROUTINE PRECIS (I,J.DELX,DELZ,NROW,NCQL,ZZERO)
c THIS SUBROUTINE PREPARES THE INPUT INTO THE
C GROUNCWATER DISPERSION MODEL
COyMCN/FINIT/NCLASS(20,20),PERKHl20t20>,PERKV(20,20),PHI(20,20),
1 JTABLE(20),KCLASS(20,2C)
COMMCN/VELOZ/VELM20,20j,.V£LV(2Q*20),.VELRt20,20) ,ANGLE 1201'20>
COKMCN/ZEROTH/C,CSTAR,CSS,II,JJ,DT,DX,DY,OCL,OCT,L,
1AI,A2,A3,A4, 61,e2,83,KFX,.CFX,,KLXf.CLX,KFY,CFYf,KLY,CLY,.
2U,VY,ZP,GP,GPC,
3XMIN,.XMAX,IDIV,YKIN»YMAX,.JDIV
COMMON/DSETS/NGWIF,NUNIF1,NLNOV1,,NUNRE1,NUNGW1,NUNRE2,NUNIF2,
1 NUNCV2,NUNOV3,NUNOV4,NUNIF3,NUNGW2fNUNGW3^NOVREl^NOVRE2iNREREl,
2 NRERE2.NOVRE3
DIMENSION C(20,20),CSTAR(20,20),CSS(20,20),U(20,20) ,VY120«0)
DIMENSION L(20,20)
DIMENSION PERMH(20,20),PERMV(20,20)
C
C SET INDEXES FOR ROWS AND COLUMNS
750
700
C
C
C.
C
850
800
JJ=NROW
SET COUNCARIES IN BOTH DIMENSIONS
XMIN»0.
XMAXsCELXMNCOL1)
YMIN=ZZ£RO
YMAX=ZZERO+OELZ»(NROW1 )
TRANSFER NODE INFORMATION TO KMATRIX
CO 700 I=*1,NCOL
CO 750 JROW=l»NROW
L( I,JROW)=»KCLASS( It JROWl
CONTINUE
CONTINUE
SET UP ARRAYS FCR VERTICAL AND HORIZONTAL
VELOCITIES
CO 300 MUNCOL '
CO 850 JROW»1,NROW
J=NRLW+lJRQW
UUtJROW)»VELHI I»J>/864CO.
VY(I>JROW)=VELVII,J)/86400.
PERMH( I, JRQWI*PERKH( J, J )
PERMVtlf JROW)»PERKV( I,J)
CONTINUE
CONTINUE
WRITE OUT DATA FOR INTERFACE WITH DISPERSION MODEL,
WRITE (6,900)
900 FORMAT ('1*,20X,•INTERFACE DATA LISTING1)
WRITE (6,910) IltJJ
910 FORMAT (215)
WRITE (6,920)
920 FORMAT <20X,'VERTICAL VELOCITIES IN FT/SEC1)
DO 925 JROW»1,NROW
WRITE 16 ,.930) ( VY( I, JROK ),, I»l,,NCOL)
930 FORMAT (10E8.2)
925 CONTINUE
WRITE (6,940)
940 FORMAT (20X,'HORIZONTAL VELOCITIES IN FT/SEC•)
CO 945 JRQW»ltNROW
WRITE (6,950) (U(I,JROW),I«1,NCOL)
950 FORMAT uoE8.2)
945 CONTINUE
WRITE (6,960)
960 FORMAT ('1«,20X,'DISPERSION NODE CLASS MATRIX*)
VEL 85
VEL 86
VEL 87
VEL 88
PREO 1
PRED
PREO
PRED
PREO
PREO
PREO 8
PR£0 9
PREO 10
PRED 11
PRED 12
PREO 13
PRED 14
PRED 15
PREO 16
PREO 17
PRED 18
PRED 19
PREO 20
PREO 21
PRED 22
PRED 23
PREO 24
PRED 25
PRED 26
PREO 27
PREO 28
PRED 29
PREO 30
PREO 31
PREO 32
PREO 33
PRED 34
PRED 35
PRED 36
PRED 37
PRED 38
PREO 39
PREO 40
PRED 41
PREO 42
PRED 43
PREO 44
PRED 45
PRED 46
PREO 47
PRED 48
PftED 49
PREO 50
PRED 51
PREO 52
PftEO 53
PREO 54
PREO 55
PRED 56
PR60 57
PREO 56
303

c
c
c
c
c
970
965
1001
CO 965 JROW=l,NROW
WRITE (6,970) ( L( I , JROW ),.!»!, NCOL )
FORMAT (2012)
CONTINUE
WRITE INTERFACE INFORMATION ON
WRITE (7..10Q1) IDIV,XMIN,XMAX
FORMAT (I5,2F10.2)
WRITE (7,1001) JDIV.YMIN.YMAX
CO 599 J=.1,JJ
WRITE (7..651) (HI, J), 1 = 1,11)
CARDS
651
599
602
601
606
605
612
611
616
615
(U(I,J), 1=1,11)
(U(
1=1, II)
(VY( I,J), 1 = 1,11)
(VY(I,J), 11,1!)
WRITE (NGWIF)
FORMAT (4012)
CONTINUE
CO 601 J=.1,JJ
WRITE (7*602)
WRITE (NGWIF)
FORMAT (10ES.2)
CCNTINUE
CO 605 J=1,JJ
WRITE (7,606)
WRITE (NGWIF)
FORMAT (10E3.2)
CONTINUE
CO 611 J=1,JJ
WRITE (7,612) (PERMHlI,J),1=1,II)
WRITE (NGWIF) (PERMh(I,J).tI»if II)
FORMAT (10E8.2)
CCNTINUE
CO 615 J=1,JJ
WHITE (7..616) IPERMV( I.J),1 = 1,II)
WRITE (iJGWIF) (PERMV( UJ),11,11)
FORMAT (10E8.2)
CCNTINUE
ENC FILE NGWIF
REWIND NGWIF
RETURN
ENC
SUEROUTINE GWQUAL
C PRGGKAM GWQUAL.
C COMPUTES CONCENTRATIONS OF A POLLUTANT AT ALL NODES
C IN A GRIC SYSTEM.
C« DEFINITION OF SYMBOLS
C* A1A4 WEIGHTING COEFFICIENTS ON THE CONVECTIVE TERMS
C« Bl83 WEIGHTING COEFFICIENTS ON THE TIME TERMS
C» CFX BOUNDARY VALUE OF CONCENTRATION
C» CFY BOUNDARY VALUE OF CONCENTRATION
C* CLX BOUNDARY VALUE OF CONCENTRATION
C» CLY BOUNDARY VALUE OF CONCENTRATION
C» GGUIT REFERENCE VALUE TO ENC CALCULATIONS
C» CCL*CIFFUSION COEFFICIENT IN THE X DIRECTION
C* DCT=DIFFUSION COEFFICIENT IN THE Y DIRECTION
C« CT=VALUE OF THE TIME STEP
C« CTICE=T1ME STEP BETWEEN DIFFERENT TIDAL VELOCITIES
C« CTMAX=MAXIMUP VALUE OF CT
C» DTMIN=MINIMUM VALUE OF CT
C« CTSTEP«MULTIPLIER FOR SUCCESSIVE TIME STEPS
C« CX=GRID SPACE INTERVAL IN THE X DIRECTION
C» CY»GRID SPACE INTERVAL IN THE Y DIRECTION
C» IDIV=.NUMBER OF GRID SPACES IN THE X DIRECTION
C» JDIV=NUMBER OF GRID SPACES IN THE Y DIRECTION
C» II=NUMBER OF GRIC POINTS IN THE X DIRECTION
C» JJ=NUMBER OF GRID POINTS IN THE Y DIRECTION
C» ICOUNT*COUNTER FOR. THE NUMBER OF TIME STEPS
C» IFREC=FRECUENCY FOR READING NEW TAPE VELOCITIES
C» ITICE*REFERENCE PARAMETER BETWEEN MAIN AND VCALC
C» JCUIT,JQUIT=SUBSCRIPTS FOR THE REFERENCE POINT
PRED 59
PRED 60
PRED 61
PRED 62
PRED 63
PRED 64
PRED 65
PRED 66
PRED 67
PRED 69A
PRED 70
PRED 71
PRED 71A
PRED 72
PRED 73
PRED 74
PRED 75
PRED 76
PRED 77
PRED 78
PRED 79
PRED 80
PRED 81
304

ISET»JSST=COURDINATES FOR POINTS OF INPUT ON BOUNDARY*
JMIN»JMAX=LIMITING POINTS FOR INPUT ON THE BOUNDARY
KIN*INTEGER SPECIFYING THE TYPE OF INPUT
SPECIFYING
SPECIFYING
SPECIFYING
SPECIFYING
KFX=INT6GER
KFY=INTEGER
KLX=INTEGER
KLY= INTEGER
NHR=INTEGER
THE
THE
THE
THE
THE
E.G.
B.C.
B.C.
B.C.
TIME
ON
ON
ON
ON
FOR
THE
THE
THE
THE
LEFT BOUNDARY
TOP BOUNDARY
RT. BOUNDARY
BOT. BOUNDARY
PRINTING
TRIDIAGONAL MATRIX
MATRIX
MATRIX
ALGORITHM
SPECIFYING
?VN=NUMBER OF TIME STEPS
BETAXtbETAY*CQEFFlCIENTS IN THE
GX,GY*C06FFICIENTS IN THE TRIDIAGONAL
ZX,ZY=COEFFICIENTS IN Tt.E TRIOIAGCJNAL
ZP,GP,GPC»VARIABLES USEC IN THE THOMAS
DEL1DEL10=R.H.S. GOEFFICIENTS
R01RC3*R.H.S. COEFFICIENTS
P1P3,0103 »SlS3*R.t.S. COEFFICIENTS
TICE=;TIME TO CALL UP A NEW TIDAL VELOCITY FIELD
TINPUT=LENGTH OF TIME OF INPUT
(j=V:LOCITY FIELC IN THE X DIRECTION
VY=veLCCITY IN THE Y DIRECTION
V£L=NAME OF THE TIDAL VELOCITY FIELD ON TAPE
LtI,J)= SPECIFIES THE GEOMETRIC CONFIGURATION
CU,J)=RESULTING COMCENTRAT ION FIELD
CSTAR(ltJ>»!NTERMEDIATe CONCENTRATION FIELD
CSSf I,J)=INTERMECIATE CONCENTRATION FIELD
C*
C»
C»
C*
C*
C*
C*
C»
C*
C*
C*
C*
C*
C*
C*
C*
c»
C*
C*
C*
c»
c»
C*
C*
C —
C
C
C
C
C
C
C
C
C
C
C
C
C —
c «****« **««««*•******»<
DIMENSION C(23,21),CSTARt23,21),CSS(23».21l»U(46.34) ,VYU6«34>
CI^ENSION L(23,21),VELU9,34)
DIMENSION GX(23J,ZX(23J
DIMENSION GY(2l I,ZY(211
DIMENSION PERMH(20,2QJ,PERMV(20,20)
DIMENSION DCLt20,20)tDCT(20,20)
DIMENSION TITLE(IO)
COMfCN/ZEROTH/C,CSTAR,CSS»II,JJ,OT,DX,DY»DCL,.DCT,L»
lAl,A2,A3f A4,Bl,a2,B3,KFXt.CFX,KLX,CLX,KFYiCFY,KLY,CLYi.
2UtVY,2P,GP,GPC,IFRSQ,ICOUNT,VELi
3XMlNrXMAXiIOIV,YMIN,YMAX,JDlV
COMMCN/FIRSr/ZX,BET4X,GXhDELl,DEL2fDEL3iDEL^,0£L5
COMMCN/SECONC/ZY,B€TAY,.GY,DEL6,DEL7tDEL8»UEL9,OEL10tROl»R02,R03
CUMMO.N/THrRO/PltP2,P3,P<»»P5,Ql,Q2«.Q3,Sl»S2»S3
COMMCN/CN/CNltCN2tCN3
COfyON/DSETS/NGWlF,NUNlFl,NUNOVl,NUNREl,NUNGHl«NUNRE2»NUNIF2»
1 NU,NOV2f KUNaV3,N'L'NOV4,NUNIF3,NUNGW2rNUNGV«3»NOVREltNUVRe2»NREREl,
2 NRERE2»NOVRE3
NUMBERING SYSTEM FOR DISPERSION NODE CLASS MATRIX KCLASS
0 = EXTERIOR
1 = LEFT BOUNDARY
2 = LEfT TOP
3 = LEFT BOTTOM

WRITE (6,2000) XMAX.DX
2000 FORMAT (' XMAX='.F10.1,5X,' DX *'»F10«AI
C GRID SIZc IN THE Y CIRECTION
READ 15,1000) JCIV,YMIN,YMAX,DY
WRITE (6r2001) YKAX.OY
2001 FORMAT (' YMAX=',.F10.4,5X, ' DY »',F10.4)
C
C
SET BOUNCARY CONDITIONS
REAC (5,2003) KFX,CFX,KLX,CLX
READ (5,2003) KFV,CFY,KLY,CLY
2003 FORMAT (2(13,F7.2))
WRITE (6,2004) KFX,CFX,KLX,CLX
2004 FORMAT!' KFX=',I3,.' CFX=',F7.2,' KLX = ',I3,« CLX=',F7.2)
WRITE (6,2035) KFY,CFY,KLY,CLY
2005 FORMAT (• KFY=',I3,« CFY»«,F7.2t' KLY*',I3,*
CONNECTIVE TERM WEIGHTING COEFFICIENTS*
REAC(5,1001) A1,A2,A3,A4
1001 FORMAT (5F10.5)
WRITE (6,2006) A1,A2,A3,A4
C
C
2006 FORMAT!' Al=l,F10.5,' A2*.«tF10.!
READ (5,103) IFREO.TMAX
103 FORMAT!I3tF10.5)
UMTS, OF TMAX ARE DAYS
A3=' ,F10.5 , «
TERMINATION OF CALCULATION
CCUIT» IQUIT.JQUIT
C CONCENTRATION FOR
REAC (5,1004)
C
C TIME STEP SECUENCE NO. OF STEPS, Mlivu MAXt STEP SIZE
REAC (5,1000) NN,DTMIN,.DTMAX,OTSTEP
WRITE (6,7008) NN,DTMIN,DTMAX,DTSTEP
2008 FORMAT I* ',15,' TIME STEPS',' OTMIN* »,F10.4,' OTMAX* f,F10.4t
1' CTSTEP= SF10.4) •
C UNITS GF OTMIN AND CTMAX MUST BE CONVERTED FROM HOURS TO SEC
CTMIN=DTMIN»360Q.
LTMAX=DTMAX*3600.
C
C
INPUT CONCENTRATION (TYPE,FROM JMIN TO JMAX)
REAC (5,1203) KIN,JMIN,JMAX,ISET^JSETtTINPUT
; UNITS CF TINPUT ARE CAYS
1200 FORMAT (5I3.F10.5)
fcRITE (6,1500) ISET,JMIN,.JMAX
1500 FORMAT (' INPUT AT I«»,I3,' BETWEEN J««,I3,' TO J»»tI3)
WRITE(6,1501) TINPUT
1501 FORMATC TINPUT* ', F10.2 )
1 SPACE INCREMENTS
DIVJaJDIV
CWI=It)IV
; NUMBER CF GRID POINTS IN X DIRECTION
II=ICIVH
; NUMBER OF GRID POINTS IN THE Y CIRECTION
JJ=JCIV+1
WRITE (6,2007) II,JJ
2007 FORMAT (« II=',I3,« JJ*',I3>
C
C
REAC AND WRITE THE L ARRAY
WRITE16.85)
85 FORMAT!/' L(I,J) ARRAY FOLLOWS')
IF (1121) 90,90,.95
90 CO 92 J=l,JJ
91 READ (5,98)!L(I,J),1=1,II)
91 REAC (NGWIF) ( L ( I ,J ), 1 = 1, II )
92 WRITE 16..99) (L
GO TO 89
95 CO 97 1=1,11
306

96 READ 15,98) (L( IRGV , J ) ,J*!, JJ >
97 WRITEl6t99) ( Ll IREV, J 1 ,.J = 1» J J )
98 FORNAT UOI2)
99 FORMAT « A« ',4012)
89 CONTINUE
C
C INITIALIZE ALL B.C. TO ZERO
CO 100 I=5li II
DO 100 J=UiJJ
CH,J)=0.
CSS( I,JJ=50.
100 CONTINUE
CSTART=0.1E30
C TEMPORARILY, USE NONZERO VALUES TO AVOID UNDERFLOWS
CO 150 Inltll
00 150 J=.1,JJ
IF (LI UJJ.EO.O) GO TO 150
C( I, J)=CSTART
CSTARtlf J)»CSTAftT
CSS( I,J)*CSTART
150 CONTINUE
C
C SET THE VELOCITY FIELD
DO 101 I'l.H
CO 101 J^ltJJ
U(I,J)=0.
VYI I, J)=0.
101 CONTINUE
C
C FIRST TIKE STEP
CT=DTHIN
C
C SET CONSTANT VALUES
110 T*0.
ICOUNT»0
TIUE=0.
CTICE600*IFREQ
CT=CT/OTSTEP
CONVERT TMAX AND TINPUT FROM DAYS TO SECONDS
TMAX=TMAX»86<»00.
TINPUT=TINPUT*86400.
B20.1667
B3=B2
C INSTEAD OF USING VCALC, REAO VELOCITY FIELD BELOW
C VELOCITIES ARE IN (FT/SEC)
DO 601 J=.ifJJ
C READ 15,602) 606) ( VYU , J > i 1 = 1, II)
READ (NGWIF) ( VYI I » J ) i I=l« II )
606 FORMAT MOES.2)
605 CONTINUE
C READ IN hORIZCNTAL AND VERTICAL PERMEABILITIES
CO 611 J=il»JJ
C REAO{5,612» ( PERMH( I » J ) , 1 = 1 , 1 1 }
READ (NGW1FMPERMH1 I,JJ,I = 1, JI)
612 FORMAT ( 10E8.21
611 CONTINUE
00 615 J«,1,JJ
C REAO (5,6161 ( PERMV( I, J J rI«UII )
READ(NGWIF) ( PERMV( I , J ) , 1*1, U )
307

616 FORMAT (1068.2)
615 CONTINUE
REHlND NGWIF
C
C COMPUTE DISPERSION COEFFICIENTS
C BASED ON HOOPES AND HARLEMAN FORMULAS
C FORMULATION BELOW COMPUTES COEFFICIENTS IN ( SQ CM/SEC)
VISKINO.C0001040
C VISKIN ABOVE IN (SQ FT/SEC), FOR T6MP 71.5 F
C CONVcRT TO (SQ CM/ScC) BELOW
VISKIN»VISKIN*0.001Q75
C BELCW CONVERT PERMEABILITY IN (GPO/SQ FT) TO ISO CM)
C ALSO* CONVERT VELOCITY FROM (FT/SEC) TO (CM/SEC)
ALPHA3»88.
BETA31.2
DO 350 I=U,II
CO 350 J^ltJJ
C USE ABSOLUTE VALUES UF VELOCITY
UABS=;ABS(U< I, J) )*30.S
VAUS = ABS(VY{ I,J) )*30.5
Ph=PeRMH(I,J)«0.542E9
PV=PERMV( I,J)»0.542E9
TERM1=VISKIN*ALPHA3
TER^2=IUABS»(PH»»0.5M/VISKIN
CCLU,J)aTERMl»(TERM2«*BETA3)
TERM3«(VABS»(PV»*0.5))/VISKIN
DCTU,J)=TERM1»(TERM3»»BETA3)
350 CONTINUE
C CONVERT CCL AND OCT FROM ( SQ CM/SEC) TO I SQ FT/SEC)
CO 360 Mlt II
DO 360 J»1,JJ
DCLU,J)=DCL< I* J)»0. 001075
CCT(I,J)=»DCT(I,J)*0.001C75
360 CONTINUE
C
C BEGIN TIPE 00 LOOP
C
400 CO 530 M4l,NN
DARG^DT
IF (CARGCTMAX) 430i4ZOi420
430 CONTINUE
C OEGIN TIKE STEPPING SEQUENCE
IF(M3) 402i402»403
402 DTf«3630.
GO TC 419
403 IF (M8) 404t405»40S
404 CT»DT«DTST6P
GO TO 419
405 IF(Mll) 406,407,408
406 CT»48.f3600.
GO TO 419
407 CT120.«3600.
GO TO 419
408 IFIM13) 404,409,420
409 DTDTMAX
GO TO 419
420 CTOTMAX
419 CONTINUE
C
C CHECK FCR SINGLE OR COUBLE INPUT
C STATEMENTS BELOW EXCLUDED TEMPORARILY
C IF (KIN. EC. 1) GO TO 433
C T143200.OT
C T243200.+TINPUT
C IF (T.LE.T2.AND.T.GE.T1) GO TO 435
433 CONTINUE
IF .(TTINPUT) 435,435,437
308

435 CONTINUE
C CALL SUB. CINPUT TO SET INPtT VALUES
CALL CINPUTUWIN,JMAX,ISET*JSET,KINfTiTINPUT)
437 CONTINUE
ICCUNTICOUNT*!
T*T+DT
C
C SAVE C(N) BY STORING IN CSS
DO 440 Itl*II
00 440 J^l.JJ
440 CSSII,J)* C( ItJ)
C
/• — — M.^•••^••B—«•* »»"•*•*• ^•••» — **"™**^—™*—<—> »»^^*^»"«
C BYPASS USAGE INSTRUCTIONS FOR SUBROUTINE VCALC
C SET UP MEW VEL FIELC AT TIDE INCREMENTS
C IF UCOUNT.EC.l) GO TO 700
C IF (T.EU.TIOE) CO TO 700
C GO TO 800
C 700 CONTINUE
C CALL VCALC (ITICE)
C CONTINUE
C 800 CONTINUE
C THREE EQUATION SCHEME FOR SOLUTION
C SOLVE THE.X EQN. FOR CSTAR
CALL TRIX
445 CCNTINUE
C RESTCRE C TO C«N) FROM CSS
CO 480 lilt II
CO 480 JiltJJ
480 C(I«J)*CSS(IrJ)
C SCLVE THE Y EQN. FOR CSS
CALL TRIY
CONTINUE
C SCLVE EXPLICITLY FOR CtN+11
CALL EXPLIC
439 CCNTINUE
C PRINT THE CONCENTRATION CIST. IF T EQUALS TIDE
C IF (ICOUNT.EQ.il GO TO 250
C IF (T.EQ.TIOE) GO TO 250
C GO TO 285
C 250 CONTINUE
C IF HITICE/NHR)«NHR.NE.IT10£) GO TO 275
C PRINT OUT RESULTS AT ALL TIKE STEPS
CALL OUTPUT!T,DTMAX)
275 CONTINUE
ITICEITICE*l
TIDETIDE+OTIDE
285 CCNTINUE
IF(TTMAX) 520r540,540
520 IF (CUQUIT,JQUIT).CQ.um 530i540f540
530 CONTINUE
540 CONTINUE
RETURN
END
SUBROUTINE VCALCtITICE)
C SUBROUTINE VCALC
DIMENSION C(23,21),CSTAR<23i21).CSS(23»2I)tUl46«34»,VY(46rf34l
DIMENSION L(23i2DtVEL(49,34)
DIMENSION DCL<20,20)»OCT<20,20)
COI"MCN/ZEROTH/C,CSTARtCSSrIItJJtDT*OXtDY»OCLtDCT»l»
lAl,A2,A3,A4,Bl,B2iB3»KFXiCFX,KLX»CLXfKFYrCFY,KLYiCLY>
2UrVY^ZP,GP,GPC,IFREQ,ICOUNTiVEL
C REAO IN VELOCITY FIELD FROM TAPE
IF

CO 18 NFREQ=1, IFREQ
18 REAC (11) VEL
IMINC = (72IS«m/iFR£Q
GO TO 45
21 CONTINUE
IF (ITIOE.EU.IW1N01 GO TO 25
DO 23 NFREQ=1, IFREQ
23 REAC (11) V6L
GO TO 45
25 REWIND 11
CO 27 NFREQ=1, IFREQ
27 REAC (11 ). VEL
IWIND=IWINDK72/IFREQ)
45 CONTINUE
18=45
i=ra
JR=33
JJM1=JR1
C TRANSFER VELOCITY FIELD FROM VEL TO U AND VY
C STORE U AND VY FOR CCDEVEN AND EVENODD POINTS
5 CO 10 J=i2tJJMl»2
U( IT J)=VEL( I, J)
VYU fJ)=VEH I, J41)
10 CONTINUE
1 = 11
IF U) 50,50,15
15 DO 20 J=U,JR,2
UU,JJ=V£L( I, J>
VY( UJ) = VEL( I,J + ir
20 CONTINUE
1 = 11
GO TC 5
50 CCNTINUE
C CALC FOR OCDODO PUINTS
I=IB
40 CO 190 J=vl,JR,2
C K CHECK FOR U AND VY
IF (K9) 55, 75^ 110
55 IF (Kl) 100.110,. 60
60 IF (K3) 120,130,65
65 IF 1K6) 145,100, 70
70 IF (K8) 100,110, 75
75 U(M,N»= .25«(Ul I1^J)+:U( I + l,4)fU( I, Jl )+U{ 1 1 J+l
VY(^^N) = .25»(VY( 11, J)4iVY(IH,J)*:'VY(
GO TC 190
100 UIM,N)»0.
VY(M^N)=0.
GO TC 190
110 U(MfN)».5*lUU,JU)+U

GO TC 190
170 U(I*,N) =
190 CONTINUE
1=12
IF (11) 350,40,40
350 CCNTINUE
CO 700 M*2,46,2
DO TOO N=i2,34i2
J = N/2
U(l,J>=U(M,N)
VY< UJ> = VYIM,N)
700 CONTINUE
WRITE CUT VELOCITY FIELD
NHR=L8/IFREQ
IF (UTICE/NHRMNHR.NE.ITIOE) GO TO 500
WKITE (6,410) TIME
VjRITL 16,420)
CO 400 IR*l» II
WRITE (6,401) I,

SM=2.
GO TO 100
70 CONTINUE
40 JM1=J1
SIGM=1.
SM = 0.
IF (JJJ) 55V55.A2
TOP CR BCTTOM
42 IF (K9) 45,100,100
45 IF (Ki) 20,100,46
46 IF (K3) 50,55,47
47 IF(K5) 50,55,48
48 IF (K7) 50,55,100
TOP CONDITION
50 JM1=J+1
GO TC (30,60,70),.KFY
GO TO 100
BOTTOM CONDITION
55 JP1=J1
GO TO (110,115,120) ,KLY
110 SIGP=1.
SP=Q.
GO TC 141
115 SIGP=1.
SP=2.
GO TO 141
120 CONTINUE
100 JPl=J+i
SP=0.
141 CONTINUE
CALL COEFX( I,J)
IF (K.EU.10) GO TO 191
IF (IM) 425,425r375
375 IF (V(I1» <»00,425,400
C LEFT OR RIGHT
400 IF (Kl) 450,425,401
401 IF (K10) 402,191,425
402 IF (K9) 403*450,191
403 IF (K^) 425,450.,404
404 IF (K6) 450,475,475
C LEFT BOUNDARY CONDITION
425 GO TO (430,431,432)tKFX
430 CIMU=C( I+1,J)
CIP1J»C( I+1,J)
GO TO 485
431 CItUJaC
CIP1J = C( I+1,J)
V(I)= BETAX*2.#ZX( I)
GO TO 485
432 CONTINUE
C R+GHT BOUNDARY CONDITION
475 GO TO (480,481, 482), KLX
480 CIP1J=C< Il.J)
CIM1J*C< Il.J)
V(I)= EETAX(ZP»GP)/V(Il)
GO TO 485
481 CIP1J=CII1, J)+2.*CIUJ)
CIM1J=»C( I1,J)
V(I)= BETAX + 2.*GX(n(ZP«GP)/V(IU
GO TC 485
*82 CONTINUE
C INTERIOR CONDITION
450 CIPlJaC( I+l.J)
312
VJI)= BETAX(2P»GP)/V(Il)
485 CONTINUE

CIJ=C(I,J)
CIJMl = SIGM«Ct It JK1.)+SK«C( I»J)
CIJP1=SIGP»C(I,JP1J+SP*C(I, J)
RHS =. SUMF(OELl,CIMiJ,OeL2,CU,DEL3,CIPlJ,DEL4,CIJMl,DEL5i
1CIJP1)
IB=l
IN=2
172 IF IIM1 186,186,187
187 IF (V(ID) 188,186,188
186 MI)=RHS
GO TO 189
188 W(II = RHStW< Il)*Zf>)/V< IU
1B9 CONTINUE
C CHECK FCR TYPE 2 B.C.
GO TO 195
191 W(1I=C(UJ)
V( I )=.BETAX
195 CONTINUE
C WRITE OUT W AND V i TEMPORARILY
C WRITE <6,.2500) I W( I ), 1 = 1, 11 )
C WRITE (6,2530) ( V( IJ , I* 1 iI I )
2500 FORMAT (10E8.2)
C BEGIN CSTAR CALCULATION
CSTARdRtJ)=W« IR)/VUR)
ISCL=IRM
00 2^0 I=?1,.ISOL
IRtV=lRI
CALL CQEFXIIREV.J)
IF (V(IKEV» 230,210,230
210 CSTAk
C SUB COEFX
DIMENSION C(23,21),CSTAR{23,21),CSSI23,2l),UU6,34) .VY146V34)
DIMENSION LI23.21)
CIMENSION DCLJ20,20),DCT(20,20)
DIMENSION GX(23)>IX<23)
CO^MCNVZcROTh/C.CSTARtCSS,!!, J J ,DT, OXtDY^CL.OCT.L,
lAl,A2,A3»A4,Bl,B2,a3,KFX,CFX,KLX,.CLX,KFY,.CFY,KLY«CUY^
2U,VY,ZP,GP^GPC,IFREQ.ICOUNT
CayMON/FlRST/ZX,aETAX,GXrDELl,DEL2,OEL3,OEL4,OEL5
CXSCDX«DX
DYSQ=.OY»CY
21 K=LII,J)
IF(K10)9,9,70
9 CONTINUE
2X(II»B3/DTDCL(I,J)/DXSO(A4»U«I,J))/OX
BETAX = B1/DT + 2.0*CCL( I,.J I /OXSQ + t A4A3)*U( I , J)/OX
GX(I)=B2/DTCCL(I,J)/OXSQ *A3*U(I«J)/OX
CEL1=B3/CT*A2*U1I,J)/DX
CEL2 = 81/DT2.0«CCT{I,J>/DYSQ*;IA1A2)#U(I,J}/DX
CEL3*Q2/CTA1*U(I,J)/OX
CEL4=DCT(I,J)/DYSQ*VY(I,J)/(2.0»DY)
10 CEL5sCCT

50 IF (L(.IltJ)lO) 51,59,57
51 IF (U I1,J> 9) 55,60,59
55 IF (L(IUJ)l) 20,57,56
56 IF (1(11,J)4) 57,60,60
57 ZP=ZX,CSS(23,21),U(46,34) ,VY146,'34)
DIMENSION L(23,21)
DIMENSION DCL(20»20),DCT(20,20)
DIMENSION GY.(21)rZYC2l)
DIMENSION V(23)iW(23)
CO^NCN/ZEROTh/C,CSTAR,CSS,IUJJ,OT,DX»DY,DCL,DCT,L,
l/Sl,A2,A3^A4,Bl,B2,B3,KFX,CFX,KLX,CLX,KFY^CFY,KLY^CLY^
2U,VY,ZP»GP,GPC,IFREQ,ICOUNT
COV^CN/SECOND/ZY,B£TAY,GY,OEL6,DEL7,OEL8,.DEL9,DELlOfR01,R02,R03
SUMF{Sl,Pl,S2,P2,S3,P3t,S4,.P4, S5, P5I »Sl
314

Kl
CD 550 I=iM, 11
J"0
20 J=J + l
If (JJJ1. 55Cu2i»Zl
21 K=L(IfJ>
IF fKU 20,140*140
140 K=J
C FIND TH£ LOWEST BOUND
5
IF (l(I, JREVl) 10,10,15
10 J=J*i
GQ TC 5
15 JBOTsJREV
DO 465 j=N,jear
IF(Jl) 320t320.3lQ
320 IWl*l*vl
SC re (330,33l»3321tKFX
330
GO TG 370
331 SIGM=,1,
SH=2.
GO TC 370
332 CONTINUE
310
IF UIl> 53,.58,42
C LEFT OR RIGHT
42 IF (K9> 45,370,370
*5 IF Ul> 20,57,
IF (K..EQ.IOJ GO TO
IF JJN) 475,475,575
575 IF (VOJ1J) 600,475,600
C TOP CR BOTTOM
600 IF 1K1) 460,460,601
60L IF

V(J)=8ETAY
GO TO 455
CUM— C(I,J + 1)*2.»CII*.J)
CIJP1»C{ I,J*1J
V( J)CETAY»2.»ZY{J)
GO TO 455
432 CONTINUE
C 8CTTCM CONDITION
485 GO TO (490,491(4921,KLY
490 CIJP1*C( ItJ1)
CIJM = Ct I.J1)
VIJ)*.SETAY tZP»GP)/VUl)
GO TO 455
491 CIJP1— C( I, Jl)+2.»Cl IiJ I.
VUMBETAY t2.»GY(J)CZP»GP)/V{ Jl)
GO TO 455
492 CONTINUE
C INTERIOR CONDITION
460 CIJP1*C( I.J+l)
CIJM1»C( I»Ji)
VUU.BETAY (ZP«GP)/V(J1)
455 CONTINUE
CIJ=C(I»J)
CIMlJaSIGM*CUMl,.J)+SM»CCLX,KFY,CFY^KLY,CLY,
2U,VY,ZP,GP,GPCt IFREQ.ICOUNT
316

COMMCN/S6COND/ZY,QETAY,GY,DEL6,DEL7,DELG»DEL9,OEL10»ft01fR02,R03
DXSC=CX»CX
DYSG»OY*OY
21 K»LCI,J)
IF (K10) 9,9,70
9 CONTINUE
ZY< J)»DCT(I,J)/DYSQA4*VY( 1, J)/DY+B3/DT
BETAY=Bl/DTf2.0«CCT(I,J>/OYSCH:.(A/DY + B2/DT
CEL6=>DCTU,J)/CYSQtU2.5>»VYU,J)/OY+B3/OT
DEL7^2.0«DCT(I,JJ/DYSQ(A2A1)«VY(I,J)/DY
DEL8*CCTU,Jl/CYSQ(Al.5)«VYt I ,J )/DY + B2/OT
CEL9=*B3/DT
CEL10—B2/DT
R01=01/DT
R02=B2/DT
10 R03*B3/DT
TOP CR BOTTOM
IF 55,.60,59
55 IFCL(IiJl)l) 20,60,56
56 IF(LII,JU3) 57,60,51
51 IFU(IrJn3) 57,60,58
58 IF(L(IrJD7) 57..60.60
57 2P=ZY{J)
GO TO (52,53,54),,KFY
52 GP=GY(J1)+ZY(J1)
GPC»GY(J)
GO TO 100
53 GP»GYUl)ZYUl>
GPOGYUl
GO TO 100
54 CONTINUE
59 ZP«ZY(J)
GP0,
GPC=GY(J)
GC TC 100
INTERIOR CONDITION
60 ZPZY(J)
GP=GY(J)
GPC»GYtJ)
GO TO 100
20 ZP»0.
BETAY0.
GP»0.
GPC'O.
GO TO 100
TOP CONDITION
34 ZP*0.
GP»0.
GO TO <35136i37>,,KFY
35 GPC=GY(J)*ZY(J1
GO TC 100
36 GPCGY(J)ZYU)
GO TO 100
37 CONTINUE
BOTTOM CONDITION
33 GP»GY(J1)
GO TO I38,39,40),KLY
38 ZP*ZY(J)^GY(J)
GPC»0.
GO TO 100
39 ZP*ZYtJ)GY(J)
317

GPC*G.
GO TC 100
40 CONTINUE
42 ZP=0.
GP=0.
GPC=0.
EETAY=1.0
GO TO 100
70 CONTINUE
100 CONTINUE
RETURN
ENC
SUbROUTINE EXPLIC
C SUBEXPLIC
DIMENSION C(23,21),CSTARt23,21),CSS<23,21),U(46,34) ,VYt46,'34)
DIMENSION L(23,21)
DIMENSION CTEMP(23)
DIMENSION DCL(20,20),DCT(20,20)
COMVCN/ZEROTH/C,CSTAR,CSS,II,JJ,DT,DX,DYrDCL,DCT,L,
1A1,A2,A3,A4,Bl,B2,B3,KFX,CFX.KLX,CLX,KFY,CFY,KUY,CLYr
2U,VY_ZP,GP,GPC,IFREG,ICCUNT
CCMtfON/THIRD/Pl,P2,P3,P4,P5,

SP=0.
GO TO 141
160 SIGP=1.
SP=2.
GO TO 141
170 CONTINUE
LOO JPiJ+1
SIOP=i.
SP=0.
141 CONTINUE
CSIGM«SIGM«C(I,JMl)+SM*C(I,J)
CSIGP*51GP«CIItJP1)+SP»C(I,JJ
CSSIGM=SIGM«CSS(I,JM1)+SM»CSS(ItJ)
CSSIGP»SIGP*CSS(IiJPl) + SP*CSSU»J)
CALL CFEXP(IiJ)
C LEFT OR RIGHT
IF IK10) 149,41, 33
149 IF (K9) 11,80.41
11 IF CK1> 241,33,12
12 IF (K4) 33,80,13
13 IF IK6) 80,34,34
C LEFT CONCITION
33 GC TC (35,37,39)»KFX
35 CTEMPUJ* P2*C< HI, J>*Ql«CSTAR+:
2S3«CSSIGP
18=1
IN=2
GO TO 85
C RIGHT CONDITION
75 CTEMP( I»*CTEMP( I)+P1»C( I,.J)*P2»C {Ii, J)*P4»CSIGW*:P5»
1 CSIGP +01»CSTAR( IltJ)+Q2*CSTAR( I»J)+S1»CSSIGM+S2»
2CSS( I,J)+S3«CSSIGP
GO TC 85
241 CTEWPID^O.
GO TO 85
41 CTEPPUUCUtJ)
85 CONTINUE
CO 210 I=«.l, II
210 CUtJI»CT6MPCl)
280 CONTINUE
RETURN
END
SUBROUTINE CFEXP(I.J)
C SUB CFEXP
DIMENSION C(23,2U,CSTAR(23,21),CSS<23,21),U(46,34),VY146,I34)
CIKENS10N LI23.2U
DIMENSION DCL(20^20"),DCT(20,20)
COfMQN/ZEROTH/C,CSTAft,CSS,II,JJ,DT,OX,DY>OCL»OCT,l.,
A2,A3f.A4, B1,B2,B3,KFX,CFX,KLX,.CLX,KFY,CFY,KLY,CLY».
319

2U,VY,ZP,GP,GPC, IFREQ, ICGUNT
COPKON/THiU.>/Pl,P2,P3,P4,P5,Ql,Q2,Q3,Sl»S2,S3
K=LU,J)
IF (K10) 75,75,100
T5 CONTINUE
OYSC=DY*OY
CTX=CT/DX
CTX2=DT/(CX*DX)
CTYsCT/DY
CTY2=DT/(CY»OY)
P1=*1.0 + U1A2)«(UU,J)»DTX*VY(I,J).*DTY)
P2*A2»U( I,J)«DTX
P3=A1»U< I,J)»DTX
P4=A2»VY< I,J)«OTY
P5=Al«VY(I,J)»RTY
G1=CCL( I,.J)»CTX2 + A4»U» I, J)*OTX
C2*2.6»CCL( I, J>»DIX2 + (A3A

GO TO 50
130 CO 140
CUSET,J).*0.
CSTAR(ISET,J)=0.
CSSl ISET»,J)=0»
140 CONTINUE
GO TO 59
INSTRUCTIONS BELOW ARE FOR TOP INPUT
150 CONTINUE
IF (TTINPUTJ 160,170,170
160 CO 165 1*1,7
JSET=8
C(1,JSET)CN1
CSTARU,JSET>»CNl
CSS(I,JSET)=CNl
165 CCIvTINUc
CO 166 I=il3rl4
CU»JSET)=*CN2
CSTAR(I,JSET>*CN2
CSS( I,JSET)«CN2
166 CONTINUE
1=15
CU,JSET)=CN3
CSrAR(I,JSET)CN3
CSS(I.JSET)=CN3
I»16
JSET=7
C(IiJSET)*CN3
CSTARdf JSET)*CN3
CSS(I,JSET)»CN3
1 = 17
JScT*.6
C(I,JS£T>.»CN3
CST4R(I,JSET)»CN3
CSS1I,JSET)»CN3
I»18
JSET6
CU,JSET)»CN3
CSTARU, JSET»»CN3
CSS(ItJScT>=CN3
I»19
JSET*8
CII,JSET)3CN3
CSTARd, JSET)»CN3
CSS(I.JSET)»CN3
GO TC 50
170 CO 175 UJMIN.JMAX
Cdf JSETJ«0.
CSTARCIf JSET)*0.
C'SS(I»JSET)»0.
175 CONTINUE
50 CONTINUE
60 RETURN
END
SUBROUTINE OUTPUTtT.DTMAX).
C SUB CUTPUT
DIMENSION C(23,21),CSTAR(23,21) iCSSU3»21 ) ,Ut 46.3A J
DIMENSION L(23,21)
DIMENSION DCL(20».20)>OCT(20»20)
COKMON/ZEROTH/C.CSTAR,CSS,II,JJ,DT^OXDY»DCLfOCT*L^
2U.VYt.ZP»GP,GPC, IFREQ.ICOUNT
C PRINT CONCENTRATIONS IN BAY
TH»T/3600.
IF (IIll) I25il25,135
125 CONTINUE
321

GC TC 620
610 WRITE(6,600) T
600 FORMAT ( 'ISSOXf 'TIME* 'fFll.At1 SECS»»«**')
GO TO 630
620 TH=T/3600.
WRITG< 6*615) T>
615 FORMAT (' l« » SOX , ' TIME SF11.4,' HRS •»«*»•)
630 CONTINUE
WSITE<6,100)
100 FORMAT (//30X,'C CONCENTRATIONS ARE1)
CO 440 J^lfJJ
440 WRITE(6,400) J , iC(11J)»1 = 1111)
400 FORMAT I • J = ', 13, 1 IE 11.3)
GO TO 2000
135 CONTINUE
GO TO 120
110 WRITE (6,300) T
GO TO 130
120 WRITE (6i200) TH
130 CONTINUE
WRITE (6,111)
111 FORMATCO',1 1*1 12 13 14
1 l»6 I«7 1=8 I»9 I»10
21' )
CO 150 J=1,JJ
150 WRITE (6*250) J, (C( I,J),l = i,U•
230 CONTINUE .
WRITE (6>112)
112 FORMATCO1,1 I»12 I»13 I»14 I»15
16 . 1=17 I«18 I»19 I«20 I»21
22')
CO 160 J=.1,JJ
160 WRITE (6,250) Jt(C(IiJ>tIa12»19)
200 FORMAT <'1',40X,' C CONCENTRATIONS AT TIME T» SF8.2,.1 HRSJ)
300 FOrtMAT Cl',40X,' C CONCENTRATIONS AT TIKE T= «,F8.2,; SECS1)
250 FORMAT C J = •, 13, HE 11. 3 )
2CTOO CONTINUE
10 CONTINUE"
RETURN
ENU
/*
//LKED.SYSIN DC *
CVERLAY ALPHA
INSERT UNZON1
OVERLAY ALPHA
INSERT UNZON2
CVERLAY ALPHA
INSERT UNZON3,UNFLO,COEF,EXTRA,VFLOW,WTA8LEtPSIVAL,PSIMID
INSERT SINK,MATRIX,TENS,CK.CSL,HUMID,GRAPH,PINE,CURVE,PPLOT
INSERT SCILQU,GROUP,MAT2,.OUTRIT,TRAN
CVERLAY ALPHA
INSERT OVLAN1
OVERLAY ALPHA
INSERT OVLAN2,AVE,WAV6tRECVrHIDURrLODUR»RQUAL,WRIT
CVtRLAY ALPHA
INSERT CVLAN3
CVERLAY ALPHA
INSERT OVLANA
CVERLAY ALPHA
INSERT RECl
CVERLAY ALPHA
INSERT REC2
CVERLAY ALPHA
INSERT REC3
CVERLAY ALPHA
INSERT GWBLCK,FHI,80UND,.TAOLE,EQUIP,VEL,PREDIS,GWQUALiVCALC
INSERT TftIX,COEFXtTRIYrCOEFY,EXPLIC,CINPUTtOtTPUT
/*
322

//GQ.FT01FC01 CC
//GO.FT02FU01 CD
//GO.FT03F001 CC
//GO.FT04F001 CC
//GC.FTObFOOl CD
//GO.FT09F001 CD
//GO.FT10FC01 CD
//GO.FTHFOOl DO
X/GO.FT12F001 CD
//GO.FT13FC01 CD
/7 SPACR=(TRKf (
//GO.FT1AF001 CD
//GO.FT15F001 CD
//GO.FT16FC01 CD
//GO.FT17F001 CC
//GO.FT18F001 CD
//GC.FT19F001 CD
//GO.FT20F001 CD
//GO.FT21F001 CD
//GC.FT22F001 CC
//GO.SYS1N CD *
SYSCA,.$PACE = (CYL 1 < 2 , 1 ) )
UNIT=SYSC/»t SPACE=( CYL t < 2, 1) )
SYSCA, SPACE = (CYL,.{ 2, 1) )
SYSCA,SPACE = (CYL» ( 2, U )
SYSCA,.SPACE = (CYL ,.< 2 1 1 ) )
UNIT=SYSOAiSPACE = (CYLtJ 2, 1) )
UNIT=SYSCA,SPACE=( CYLr< 2 t L»)
UNIr = SYSCAfSPACEa(CYL^^2t D)
UNI T = SYSDA, SPACE*tCYL t < 2 . 1 ) )
DSN= ALAN. APQPKA4,UNIT = 23U,DI SP"< ,CATLG) ,
20i5),RLSE)tVQL=SER»USER01,OCB*(RECFM*VS,BLKSIZE»920)
Ui\l T = SYSCAf SPACE = (CYL t < 2, 1 ) )
UNI T=SYSCA,SPACE = (CYU< 2, 1 ))
UN IT=SYSCA,SPACE= t CYL t ( 2, 1) )
UNIT = SYSCA, SPAC£ = ( CYLr( 2, 1) )
UN.I T=SYSCA, $PACE = (CYLt< 2, 1 ) )
UNIT=SYSDA,SPACE=tCYL»< 2 , 1) )
SYSCA,SPACE=( CYL , \ 2, L) )
SYSCA,SPACE = «CYL»< 2t I) )
DSN*ALAN. APOPKAA,DISP=OLD»,UNI T*2314» VOL»SERUSER01
OO.S. GOVERNMENT PRINTING OFFICE: »74 582413/76 13
323

1
5
Accession Number
2
Organization
Department of
Subject field it Group
SELECTED WATER RESOURCES ABSTRACTS
INPUT TRANSACTION FORM
Environmental Engineering Sciences
University of Florida
Gainesville. Florida 32611
Title
A Water Quality Model for a Conjunctive SurfaceGroundwater System
10
Authors)
Perez, Armando I.
Huber, Wayne C.
Heaney, James P.
Pyatt, Edwin E.
ix. Protect Destination
2] Note
EPA 16110 GEW
22
Citation
Environmental Protection Agency report number, EPA600/57U013,
May
23
Descriptors (Starred First)
*Surfacegroundwater relationships, *Surface runoff, *Groundwater flow,
*Unsaturated flow, *Water pollution, Path of pollutants, Agricultural watershed,
Florida, Nutrients, Nitrates, Phosphates, Hydrologic models, Mathematical models
25
Identifier* (Started Pint)
'Agricultural pollution, Lake Apopka Florida, Water quality models
27
increasing pressures toward efficient decisionmaking in water pollution con
trol are creating the need for improved pollutant routing models, Such mathematical
models help understand cause and effect relationships between sources of pollutants
and their ensuing concentrations at various locations in a basin.
Considered in this study were both flow and water quality processes occurring
on the ground surface, in the unsaturated soil zone and in the saturated or ground
water zone. The objective was to improve already available formulations for the
above processes and subsequently to develop a methodology for interfacing the
individual models.
Emphasis was placed on the modeling of agricultural pollution. For this rea
son, nitrogen and phosphorous were the main substances considered. The selection of
the Lake Apopka basin in Central Florida as the study area wds made in accordance
with these project goals. The data base for this basin was limited, but better
than that for other candidate areas in Florida. However, this data base appeared
typical of that available to engineers in the analysis of a practical problem. The
scope of this particular study precluded additional sampling. Current data limita
tions precluded a complete verification of the model, However, various general
conclusions could be drawn. In the future, the formulation could be viewed as an
instrument for structuring data gathering efforts.
Abstractor
Wayne C. Huber
Institution
University of Florida
WRIIC
