MEASUREMENT AND SIMULATION OF SOIL WATER STATUS
UNDER FIELD CONDITIONS
By
KENNETH COY STONE
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
V
1937
ACKNOWLEDGEMENTS
•The author would like to express his appreciation to Dr. Allen G.
Smajstrla for his guidance, assistance and encouragement throughout th
research. The author would like to thank his other committee members
for their support and guidance throughout this research and in
coursework taken under their supervision.
The author also expresses his appreciation to the Agricultural
Engineering Department for the use of its research and computing
facilities. Appreciation is also expressed to other faculty and staff
members in the department for their assistance.
Special gratitude is expressed to the authors family for their
continuous support and encouragement throughout his studies.
Finally, the author would like to express his very special
graditude to his wife, Carol, for her drafting expertise, patience and
continuous support and encouragement.
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS Ü
LIST OF TABLES v
LIST OF FIGURES vi
ABSTRACT x
CHAPTERS
I INTRODUCTION 1
II REVIEW OF LITERATURE 4
Measurement of Soil Water Status 4
Soil Water Content Measurement 4
Soil Wafer Potential Measurement 8
Automation of
Soil Water Potential Measurements 11
Water Movement in Soils 15
Soil Water Extraction 20
III METHODS AND MATERIALS 26
Equipment 26
Instrumentation Development 27
Field Data Collection 30
Model Development 33
One-Dimensional Model 36
Two-Dimensional Model 41
Soil Water Extraction 50
IV RESULTS AND DISCUSSION 54
Instrumentation Performance 54
Field Data Collection 60
Model Verification 70
Model Operation with Field Data 94
One-Dimensional Model 102
Two-Dimensional Model 109
Model Application 116
V SUMMARY AND CONCLUSIONS 124
Microcomputer-based Data Acquisition System 124
Field Data 124
Modeling Soil Water Movement and Extraction 125
Model Applications 125
APPENDICES
A. Alternating Direction Implicit Finite Differencing. . 126
B. Listing of Soil Water Movement and Extraction Model. .137
REFERENCES 163
BIOGRAPHICAL SKETCH 167
LIST OF TABLES
Table 3-1. Microcomputer-based data acquisition system
components and approximate costs 28
Table 4-1. One-dimensional distribution of water extraction
for a 15 day drying cycle for young citrus trees
with and without grass cover 68
Table 4-2. Two-dimensional distribution of water extraction
for a 15 day drying cycle for young citrus trees
with grass cover 69
Table 4-3. Two-dimensional distribution of water extraction
for a 4 day drying cycle for young citrus trees
with and without grass cover 73
Table 4-4. One-dimensional distribution of water extraction
for a 4 day drying cycle for a grass cover crop
at water depletion levels of 20 kPa and 40 kPa . ... 76
v
LIST OF FIGURES
Figure 2-1. A typical soil water characteristic curve 5
Figure 3-1. Layout of lysimeter system 31
Figure 3-2. Details of individual lysimeter soil water status
monitoring system 32
Figure 3-3. Location of tensiometers in 1985 field experiment
for young citrus trees in grassed and bare soil
lysimeters 34
Figure 3-4. Location of tensiometers in 1986 field experiment
for young citrus trees in grassed and bare soil
lysimeters 35
Figure 3-5. Schematic diagram of the finite-difference grid
system for the one-dimensional model of water
movement and extraction 38
Figure 3-6. Schematic diagram of the finite-difference grid
system for the two-dimensional model of water
movement and extraction 44
Figure 3-7. Effect of the relative available soil water
and potential soil water extraction rate on the
soil water extraction rate 53
Figure 4-1. Calibration curve of output voltage versus
pressure applied for pressure transducer No. 1. . . 55
Figure 4-2. Calibration curve of digital units versus
voltage applied for the analog-to-digital
circuit used 57
Figure 4-3. Comparisons of mercury manometer manually-read
and pressure transducer automatically-read
tensiometer water potentials during drying cycles
for 2 tensiometers in the laboratory 58
Figure 4-4. Comparisons of mercury manometer manually-read and
pressure transducer automatically-read tensiometer
water potentials for all laboratory data 59
Figure 4-5. Evaluation of pressure transducer-tensiometer No. 1
by comparison of mercury manometer manually-read
and pressure transducer automatically-read water
potentials in the field 61
vi
Figure 4-6. Evaluation of pressure transducer-tensiometer No. 2
by comparison of mercury manometer manually-read
and pressure transducer automatically-read water
potentials in the field 62
Figure 4-7. Evaluation of pressure transducer-tensiometer No. 3
by comparison of mercury manometer manually-read
and pressure transducer automatically-read water
potentials in the field 63
Figure 4-8. Evapotranspiration rate from the 1985 field
experiment with a 20 kPa soil water potential
treatment young citrus tree with grass cover. ... 66
Figure 4-9. Evapotranspiration rate from the 1985 field
experiment with a 20 kPa soil water potential
treatment young citrus tree with bare soil 67
Figure 4-10. Evapotranspiration rate from the 1986 field
experiment with a 20 kPa soil water potential
treatment young citrus tree with grass cover. ... 71
Figure 4-11. Evapotranspiration rate from the 1986 field
experiment with a 20 kPa soil water potential
treatment young citrus tree with bare soil 72
Figure 4-12. Evapotranspiration rate from the 1986 field
experiment with a 20 kPa soil water potential
treatment with a grass cover 74
Figure 4-13. Evapotranspiration rate from the 1986 field
experiment with a 40 kPa soil water potential
treatment with a grass cover 75
Figure 4-14. Soil water potential-soil water content
relationship for Rehovot sand 79
Figure 4-15. Hydraulic conductivity-soil water content
relationship for Rehovot sand 80
Figure 4-16. Simulated results of soil water content profiles
for infiltration into a Rehovot sand under
constant rain intensity of 12.7 mm/hr 81
Figure 4-17. Simulated results of soil water content profiles
for infiltration into a Rehovot sand under
constant rain intensity of 47 mm/hr 82
Figure 4-18. Soil water potential-soil water content
relationship for Yolo light clay 83
Figure 4-19. Hydraulic conductivity-soil water content
relationship for Yolo light clay 84
Figure 4-20. Soil water potential-soil water content
relationship for Adelanto loam 85
Figure 4-21. Hydraulic conductivity-soil water content
relationship for Adelanto loam 86
Figure 4-22. Simulated results of soil water content profiles
for infiltration into a Yolo light clay with
initial pressure potential at -66 kPa 87
Figure 4-23. Simulated results of soil water content profiles
for infiltration into a Yolo light clay with
initial pressure potential at -200 kPa 88
Figure 4-24. Simulated results of soil water content profiles
for infiltration into a Adelanto loam with
initial pressure potential at -66 kPa 89
Figure 4-25. Soil water potential-soil water content
relationship for Nahal Sinai sand 90
Figure 4-26. Hydraulic conductivity-soil water content
relationship for Nahal Sinai sand 91
Figure 4-27. Simulated results of dimensionless soil water
content distribution for infiltration into a
Nahal Sinai sand at 57 min of infiltration time. . 92
Figure 4-28. Simulated results of dimensionless soil water
content distribution for infiltration into a
Nahal Sinai sand at 297 min of infiltration time. . 93
Figure 4-29. Two-dimensional model simulated results of soil
water content profiles for infiltration into a
Rehovot sand under constant rain intensity
of 12.7 mm/hr 95
Figure 4-30. Two-dimensional model simulated results of soil
water content profiles for infiltration into a
Rehovot sand under constant rain intensity
of 47 mm/hr 96
Figure 4-31. Evapotranspiration rate from the second drying
cycle of the 1985 field experiment with a 20 kPa
soil water potential treatment young citrus tree
with grass cover 98
Figure 4-32. Distribution of evapotranspiration through the
daylight hours 99
Figure 4-33. Soil water potential-soil water content
relationship for Arredondo fine sand. .
vi ii
100
Figure 4-34. Hydraulic conductivity-soil water potential
relationship for Arredondo fine sand 101
Figure 4-35. Simulation results for the one-dimensional model
with field data for the 150 mm depth 104
Figure 4-36. Simulation results for the one-dimensional model
with field data for the 300 mm depth 105
Figure 4-37. Simulation results for the one-dimensional model
with field data for the 450 mm depth 106
Figure 4-38. Simulation results for the one-dimensional model
with field data for the 600 mm depth 107
Figure 4-39. Simulation results for the one-dimensional model
with field data for the 900 mm depth 108
Figure 4-40. Simulation results for the two-dimensional model
with field data for the 150 mm depth Ill
Figure 4-41. Simulation results for the two-dimensional model
with field data for the 300 mm depth 112
Figure 4-42. Simulation results for the two-dimensional model
with field data for the 450 mm depth 113
Figure 4-43. Simulation results for the two-dimensional model
with field data for the 600 mm depth 114
Figure 4-44. Simulation results for the two-dimensional model
with field data for the 900 mm depth 115
Figure 4-45. Soil water potentials for the three irrigation
treatments at the 300 mm depth 118
Figure 4-46. Soil water potentials for the three irrigation
treatments at the 900 mm depth 119
Figure 4-47. Soil water storage for the three irrigation
treatments 120
Figure 4-48. Cumulative irrigation for the three irrigation
treatments 121
Figure 4-49. Cumulative drainage from the soil profile for
the three irrigation treatments 122
IX
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
MEASUREMENT AND SIMULATION OF SOIL WATER STATUS
UNDER FIELD CONDITIONS
By
KENNETH COY STONE
May 1987
Chairman: Allen G. Smajstrla
Major Department: Agricultural Engineering
Irrigation of agricultural crops is one of the major uses of fresh
water in Florida. Supplemental irrigation is required in Florida
because the predominant sandy soils have very low water holding
capacities, up to 70» of the annual rainfall occurs during months when
many crops are not grown, and rainfall is not uniformly distributed.
Under these conditions the scheduling of irrigation becomes extremely
important. The timely application of irrigation water can result in
increased yields and greater profits, while untimely applications could
result in decreased yields and profits, and leaching of nutrients.
A low cost microcomputer-based data acquisition system for
continuous soil water potential measurements was developed. The system
consisted of tensiometer-mounted pressure transducers, an analog-to-
digital converter and a portable microcomputer. The data acquisition
system was evaluated under laboratory and field conditions. Excellent
x
agreement was obtained between soil water potentials read with the data
acquisition system and those read manually using mercury manometers.
Experiments were conducted utilizing the data acquisition system to
monitor soil water potentials of young citrus trees in field lysimeters
to determine model parameters and to evaluate model performance.
Evapotranspiration, soil water extraction rates and relative soil water
extractions were calculated from the soil water potential measurements.
Two numerical models were developed to study the movement and
extraction of soil water. A one-dimensional model was developed to
describe a soil profile which is uniformly irrigated such as one which
is sprinkler irrigated. The model was used to simulate soil water
movement and extraction from a young citrus tree with grass cover in a
field lysimeter. Simulated results were in excellent agreement with
field observations.
A two-dimensional model was developed to describe a nonuniformly
irrigated soil profile as in trickle irrigation. The model was used to
simulate the two-dimensional movement and extraction of soil water from
a young citrus tree with grass cover in a field lysimeter. The model
results were in excellent agreement with field observations.
xi
CHAPTER I
INTRODUCTION
Irrigation of agricultural crops is one of the major uses of fresh
water in Florida. It was reported (Harrison et al., 1983) that in 1981,
41 percent of the fresh water use in Florida was for irrigation of more
than 800,000 hectares of agricultural crops. Three reasons cited for
supplemental irrigation in Florida were 1) sandy soils have very low
water holding capacities, 2) up to 70 percent of the annual rainfall
occurs during months when many crops are not grown, and 3) rainfall is
not uniformly distributed even during months of high rainfall.
Because no control can be placed on when and where rainfall occurs,
researchers must focus their attention on managing the soil water
content in the plant root zone. One method of managing the water
content in soils is to apply small amounts of irrigation water at
frequent intervals. In the absence of rainfall, this provides
agricultural plants with adequate water for growth and also minimizes
losses from percolation below the root zone. With rainfall
interactions, irrigation scheduling is more difficult; the objective is
to minimize irrigation inputs (maximize effective rainfall) while
optimizing production returns. Because of the complexity of these
interactions, numerical models are useful tools to study them. Many
researchers have developed models to use in the study of soil water
management (Smajstrla, 1982; Zur and Jones, 1981). These models allow
many different irrigation strategies to be investigated without the cost
normally associated with field experiments. These models involve
1
2
analyses of infiltration, redistribution, evaporation and other factors
which affect the movement and uptake of soil water by plants.
The objective of this research was to develop a data collection
system and a numerical model which together will provide researchers
with information needed for developing and validating models of soil
water movement, crop water use and evapotranspiration. The
instrumentation system developed will record soil water potentials on a
real time basis in order to provide input data for model verification
and validation. A microcomputer based instrumentation system will allow
the computer to monitor inputs and make decisions based on the input
data. The microcomputer will have the ability to monitor and control
events in the field.
Two numerical models will be developed to study the movement and
extraction of soil water. A one-dimensional model will be developed to
describe the movement and extraction of soil water in a soil profile
which is uniformly irrigated. A two-dimensional model will be developed
to describe the movement and extraction of soil water from a soil
profile which is not uniformly irrigated such as trickle irrigation.
The specific objectives of this research were:
1. To develop and test instrumentation for the real-time
monitoring of soil water potential under agricultural crops.
2. To record the dynamics of soil water movement and water
extraction under irrigated agricultural crops.
3. To develop a numerical model to simulate the soil water
extraction patterns observed under agricultural crops as a
function of irrigation management practices and climatic
demands.
3
4. To demonstrate the use of the numerical model in evaluating
and recommending irrigation strategies.
CHAPTER II
REVIEW OF LITERATURE
Measurement of Soil Water Status
The measurement of soil water can be classified into two
catagories: (1) the amount of water held in a given amount of soil (soil
water content), and (2) the potential, or tension with which the water
is held by the soil (soil water potential). These properties are
related to each other (Figure 2-1) and describe the ability of a soil to
hold water available for plant growth.
Soil Water Content Measurement
Several methods are available to measure soil water content. The
gravimetric method is the standard method of determining the soil water
content. This method consist of physically collecting a soil sample
from the field, weighing it, and oven drying the sample to constant
weight at 105 C. The difference in weights before and after drying is
the amount of water removed from the sample. Water contents can be
calculated on a weight basis, or on a volume basis if the soil volume or
bulk density was measured when the sample was taken. An advantage of
the gravimetric method is that it requires no specialized equipment.
Disadvantages are that it is laborious and time consuming. Destructive
sampling is required and sampling may disturb a location sufficiently to
distort results. This method does not lend itself to automation.
Neutron scattering allows the nondestructive measurement of soil
water content. A neutron moisture meter may be used in the field to
4
SOIL WATER POTENTIAL (mm)
Figure 2-1.
SOIL WATER CONTENT (mm/mm)
A typical soil water potential-soil water content curve.
6
rapidly and repeatably measure water contents in the same location and
depth of soil. An access tube must be installed in the soil to allow a
probe to be lowered to the desired soil depth for measurements.
The neutron scattering method operates on the principle of nuclear
thermalization. Fast neutrons are emitted from a radiation source
located on a probe which is lowered into the soil. These neutrons lose
energy as they collide with hydrogen atoms in the soil and are slowed or
thermalized. Thermalized neutrons are counted by a detector which is
also located on the probe. The number of slow neutrons is an indirect
measure of the quantity of water in the soil because of the hydrogen
atoms present in water. Because the neutron method counts only hydrogen
atoms, the method must be calibrated for each specific soil type and
location. This is especially true for soils with variable quantities of
organic (natter because the neutron meter will also count hydrogen atoms
in the organic matter and cause errors in the water content
measurements.
With the neutron method, a spherical volume of soil is sampled as
neutrons are emitted from the source. The radius of the sphere varies
with the soil water content from a small radius for wet soils to 20 or
30 cm for very dry soils. This measurement over a relatively large
volume is an advantage for homogenous soils with no discontinuities.
However when discontinuities such as water tables exist, their exact
locations cannot be detected accurately. Aribi et al. (1985)
investigated the accuracies of neutron meters when measuring water
contents near boundaries. They found that significant errors occurred
when readings were taken within the top 0.4 meters in an unsaturated
soil profile.
7
Another disadvantage of the neutron method is that the equipment
cost is high. The radioactive material contained in the neutron meter
is also a disadvantage because of the health hazards associated with
radiation. The automation of a neutron meter would be very expensive
and extremely complicated.
Gamma ray attenuation is a method which can be used to determine
the soil water content. This method determines soil water content by
measuring the amount of gamma radiation energy lost as a radiation beam
is directed through a soil. The method depends on the fact that gamma
rays lose part of their energy upon striking another substance, in this
case soil water. As the water content changes, the amount of
attenuation will also change.
Various types of gamma radiation instruments are available. One
type is intended for laboratory or stationary use. This instrument has
the source and detector parallel in line at a fixed distance apart, and
a soil column is placed between them. The source and detector are then
moved along the length of the soil column to determine the water content
distribution in the column. This would not be practical for field
application.
A second type of gamma radiation instrument has been developed for
field application. This instrument requires that two parallel access
tubes be installed into the soil. The source is placed in one access
tube and the detector in the other. Radiation is focused into a narrow
beam between the source and detector, and soil water content can be
measured in very thin layers of soil.
Another type of gamma radiation instrument consists of a stationary
detector located at the soil surface and a source located on a probe
8
which is lowered into the soil. The relationship between the source and
detector is known, and thus changes in density with depth can be
measured accurately with this instrument. Because the detector is
located at the soil surface, the instrument is limited to the top 30 cm
of the profile.
Gamma radiation instruments are expensive, and they can be
hazardous due to the radioactive source. This type of instrument would
be very difficult and expensive to automate.
Additional research has been conducted to utilize other soil
properties which would lend themselves to rapid methods of soil water
content determination. Fletcher (1939) conducted work on a dielectric
method of estimating soil water content. He used a resonance method to
determine the dielectric constant of the soil. This type of instrument
consists of an ocillator circuit and a tuned receiving circuit. The
instrument is placed in the soil and allowed to equilibrate. A variable
capacitor in the receiving circuit is then tuned to resonance, and the
capacitor reading is correlated with the soil water content. This
instrument must be calibrated for each soil in which it is used. Such a
device is not known to be commercially available.
Most of the methods discussed produce indirect measurements of soil
water content. Soil water content is calibrated to other factors such
as neutron thermalization, gamma ray attenuation or dielectric
properties. Only the gravametric method yields a direct measure of soil
water content.
Soil Water Potential Measurement
An alternative to the measurement of soil water content is the
measurement of soil water potential. A tensiometer measures the
9
potential or tension of water in the soil. The tensiometer consists of
a closed tube with a ceramic cup on the end which is inserted into the
soil and a vacuum gage or manometer to measure the water potential in
the tensiometer tube. The tube is filled with water, closed and allowed
to equilibrate with the soil water potential. As the soil dries, water
in the tensiometer is pulled through the ceramic cup. The soil water
potential which pulls water through the ceramic cup is registered on the
vacuum gage or manometer. This force is also the hydraulic potential
that a plant would need to exert to extract water from the soil.
Therefore, a tensiometer measures the energy status of water in the
soil. Tensiometers left in the soil for a long period of time tend to
follow the changes in the soil water potential.
The hydraulic resistance of the ceramic cup, the surrounding soil,
and the contact between the cup and soil cause tensiometer readings to
lag behind the actual tension changes in the soil. Lags are also caused
by the volume of water needed to be moved through the cup to register on
the measuring device. The useful range of tensiometers is from 0 to -80
kPa. Below -80 kPa air enters through the ceramic cup or the water
column in the tensiometer breaks, causing the tensiometer to fail. This
measurement limitation is not serious for irrigated crops on sandy soils
because most of the available water for plant use occurs between 0 and
-80 kPa.
Another method of measuring soil water potential is the thermal
conductance method. The rate of heat dissipation in a porous material
of low heat conductivity is sensitive to the water content in the porous
material. When in contact with a soil, the water potential in the
porous material tends to equilibrate with the soil water potential.
10
Phene et al. (1971) developed an instrument to measure soil water
potential by sensing heat dissipaton within a porous ceramic. The water
potential of the porous ceramic was measured by applying heat at a point
centered within the ceramic and measuring the temperature rise at that
point. Soil water potential measurements were obtained by taking two
temperature readings. The first one was taken before the heating cycle
and the second after the heating cycle. The difference between the two
temperature measurements was the change in temperature at the center of
the sensor due to the heat applied. The magnitude of the difference
varies depending on the water content of the porous block. Phene stated
that the sensor should measure the soil water potential regardless of
the soil in which it is embedded.
In experimental applications of the heat pulse device, the accuracy
of the sensor was - 20 kPa over the range 0 to -200 kPa. In Florida
sandy soils such a wide range of variability would not be acceptable
because most of the available soil water is contained in the soil
between field capacity and -20 kPa. Calibration for individual soils
couTd improve the accuracy of the instrument. This instrument has the
capability of allowing automated data collection and is nondestructive.
The cost of the individual sensors and the associated data recording
devices is relatively high.
Soil psychrometers are instruments which can be used to measure
soil water potentials in the range of -1 to -15 bars. They operate by
cooling a thermocouple junction which is in equilibrium with the soil to
the point of water condensation on the junction, and then measuring the
junction temperature as water is allowed to evaporate. Thus the
temperature depression due to evaporation can be related to soil water
11
tension. Only small responses are obtained at potentials above -1 bar,
making this instrument unsuitable for irrigation scheduling on sandy
soils where irrigations would typically be scheduled at much greater
potentials.
Electrical resistance blocks can be correlated with soil water
potential. The blocks are placed in the soil and allowed to
equilibrate with the soil water. Resistance blocks usually contain a
pair of electrodes embedded in gypsum, nylon or fiberglass.
Measurements with resistance blocks are sensitive to the electrolytic
solutes in the fluid between the electrodes. Thus resistance blocks are
sensitive to variations in salinity of soil water and to temperature
changes. Temperatures must be measured and resistance readings
appropriately corrected.
Resistance blocks are not uniformly sensitive over the entire range
of soil water content. They are more accurate at low water contents
than at water contents near field capacity. Because of this limitation,
resistance blocks can be used to complement tensiometers to measure
soil water potentials below the -80 kPa range. They are relatively
inexpensive, and an automated data collection system may be built around
these instruments for measuring water potentials in drier soils.
Automation of Soil Water Potential Measurements
The instruments discussed may all be used to measure the status of
soil water. Most of the techniques relate a soil property to another
property which is measured by the instrument. Two of the methods allow
a direct measure of a soil property. The gravimeteric method gives a
direct measure of soil water content and tensiometers give direct
measurement of soil water potential. The gravimetric measurement method
12
requires destructive sampling and does not lend itself to automation.
Tensiometers may be automated by recording the changes in soil water
potential, and they do not require destructive sampling. Therefore, to
measure the status of soil water under a growing crop, an automated data
collection system using tensiometers to measure soil water potential was
chosen.
To determine the soil water potential, the water potential within
the tensiometer is measured. The tensiometer water potential is assumed
to be in equilibrium with the soil water potential. Early methods of
measuring tensiometer water potentials used mercury manometers. Later,
mechanical vacuum gauges were used. Both methods functioned well for
manual applications. Neither, however, was readily automated.
Recent interest in better understanding the dynamics of soil water
movement, and in the development of numerical models to simulate this
process, has resulted in the need for an automated system of recording
tensiometer readings on a continuous real time basis.
Van Bavel et al. (1968) used a camera to take periodic photographs
of a tensiometer manometer board. By analyzing the photographs they
were able to record changes in potentials. This procedure was, however,
laborious and did not provide continuous soil water potential records.
Enfield and Gillaspy (1980) developed a transducer which measured
the level of mercury in a mercury manometer tensiometer. The principle
of operation of their transducer was the same as a concentric capacitor.
The level of mercury in the manometer corresponds directly to the length
of a capacitor plate. A steel tube was used as the outer capacitor
plate around a column of mercury. The nylon tube which contained the
mercury column acted as the dielectric material. The capacitance of the
13
transducer was measured and converted to length of the column of
mercury. This instrument was found to be very sensitive to temperature
fluctuations. Further research is needed on this transducer to make it
suitable for applications in automated data collection systems.
Fitzsimmons and Young (1972) used a tensiometer-pressure transducer
system to study infiltration. Their system consisted of a pressure
transducer which was connected to many tensiometers by a system of fluid
switches. This system required substantial time lags after switching a
tensiometer to the transducer before an accurate reading could be made.
This was primarily due to changes in volume due to the elasticity of the
system. This resulted in a large scanning time in order to read all
tensiometers. Long and Hulk (1980) used a similar system.
Bottcher and Miller (1982) developed an automatic manometer
scanning device which read and recorded mercury levels in manometer-type
tensiometers. The system consisted of a computer-control 1ed chain-drive
mechanism which moved photocells and light sources up and down the
manometer tubes. When the scanner passed a mercury-water interface, a
change in voltage was detected by the computer. This system is
expensive, uses specially constructed rather than commercially-available
components, and is not readily expandable. It also requires the use of
a central manometer board with connecting tubing from the various
tensiometers. It has the advantage of easy verification of computer
readings by manual reading of the manometers.
Marthaler et al. (1983) used a pressure transducer to read
individual tensiometers. The upper end of the tensiometers was closed
off with septum stoppers to provide air-tight seals through which a
needle connected to a pressure transducer was inserted. The needle was
14
inserted into a pocket of entrapped air in the tensiometer, and the
pressure transducer output was obtained immediately. This method
introduced an error because the air pocket was compressible and was
affected by the addition of air at atmospheric pressure in the needle,
transducer, and connectors. This system was designed to minimize errors
by minimizing the volume of air introduced into the tensiometer. The
tensiometer air pockets were, however, temperature sensitive and
introduced diurnal time lags because of expansion and contraction due to
diurnal temperature changes. Also, to read several instruments, the
pressure transducer was required to be moved manually. Thus, this type
of system does not lend itself to automation.
Thomson et al. (1982) used individual pressure transducers on
tensiometers with all air purged from the system to monitor soil water
potentials. Thus, they were able to avoid lag times associated with air
pockets in the tensiometers. They were able to read pressure
transducers very rapidly by electronic rather than hydraulic switching.
A soil water potential monitoring system that used pressure
transducers as used by Thomson et al. (1982) would be able to record
data from a number of sensors very rapidly. A continuous data
acquisition system based on tensiometer-mounted pressure transducers
would provide data necessary for models of movement of water in soils,
crop water use, and evapotranspiration. The use of a microcomputer to
monitor the pressure transducers would allow the system to be programmed
to make decisions based on the input data (Zazueta et al., 1984).
Also, the cost of a dedicated data acquisiton system would be greater
than that of a microcomputer-based system because of current
microcomputer costs and availability.
15
Water Movement in Soils
The general equations governing unsaturated flow in porous media
are the continuity equation and Darcy's Law. Hi11 el (1980) presented a
combined flow equation which incorporates the continuity equation
36
_= - v-q - S (1)
3t
where q = flux density of water,
0 = volumetric water content,
v = differential operator,
t = time, and
S = a sink or source term,
with Darcy's equation for unsaturated flow
q = - K(h) vH (2)
where K = hydraulic conductivity,
H = the total hydraulic head and defined as
H = h - z, and
h = capillary pressure head.
The resulting combined flow equation for both steady and transient
flows is also known as Richards equation
39
— = v*(K(h>v H ) - S (3)
3t
For one-dimensional vertical flow, equation (3) becomes
36 3 3h 3K(h)
— = — c K(h) — ) - - S(z,t) (4)
3t 3Z 3Z 3z
where z = the vertical dimension.
16
The sink term is used to represent the loss or gain of water from
the soil by root extraction or by application of irrigation water from a
point source.
For implicit solutions, equation (4) must be written in terms of
only one variable, soil water content or potential. By introducing the
specific water capacity, C, defined as
d e
C = (5)
d h
and using the chain rule of calculus, equation (4) may be written in
terms of the soil water potential as
9h 3 3 h 3 K(h)
C = (K(h) ) - - S(z,t) (6)
3 t 3 Z 3 Z 3 Z
For two-dimensional flow, equation (3) becomes
3 h 3 ah 3 3 h 3 K(h)
C — = —(K(h)—) + —(K(h)—) - S(x,z,t) (7)
at 3 X 3 X 3 Z 3 Z 3 Z
Th—4
o_
LO
ZZ
c
ce
i—
o
a.
c
6 -
235
~1—
240
245
250
JULIAN DAY
Figure 4-8. Evapotranspiration rate from the 1985 field experiments with a
20 kPa soil water potential treatment young citrus tree with
grass cover.
cn
cn
Figure 4-9. Evapotranspiration rate from the 1985 field experiment with
a 20 kPa soil water potential treatment young citrus tree
with bare soil.
68
Table 4-1. One-dimensional distribution of water extraction for a 15 day
drying cycle for Young Citrus trees with and without grass
cover.
Depth
Depth
Water
Water
Relative
Increment
Extraction
Extraction
Water
Extraction
(ITTT1)
(mm)
(mm)
(rnn/mm)
(%)
Tree
with Grass Cover
0-375
375
9.8
0.026
32.3
375-750
375
7.3
0.019
23.8
750-1050
300
8.7
0.029
28.6
below 1050
4.7
0.016
15.2
Tree with
No Grass Cover
0-375
375
5.9
0.016
27.9
375-750
375
7.4
0.020
35.3
750-1050
300
7.0
0.024
33.5
below 1050
0.7
0.0023
3.3
69
Table 4-2. Two-dimensional distribution of water extraction for a 15
day drying cycle for young Citrus with grass cover.
Depth Depth
Increment
(mm) (mm)
Radi us
(rrm)
100
300
500
700
Relative Water Extraction
•
( % )
0-375
375
8.59
10.46
6.92
7.66
375-750
375
5.65
5.79
5.46
5.65
750-1050
300
8.03
8.26
3.94
4.07
below 1050
4.21
3.45
6.33
6.04
Water Extraction
(mm)
0-375
375
2.6
3.2
2.1
2.3
375-750
375
1.7
1.8
1.7
1.7
750-1050
300
2.4
2.5
1.2
1.2
below 1050
1.3
1.1
1.9
1.8
Water Extraction
(mm/mm)
0-375
375
0.0070
0.0085
0.0056
0.0062
375-750
375
0.0046
0.0047
0.0044
0.0046
750-1050
300
0.0082
0.0084
0.0040
0.0041
below 1050
0.0043
0.0035
0.0064
0.0061
24. Differences between the two treatments showed that ET of the
tree with grass cover was greater than that of the tree with no cover,
especially early in the drying cycle. The 1220 mm depth tensiometer for
the tree with no grass did not function properly so these data were not
included in the data anaylysis, but this did not introduce a large
error, from Table 4-1.
The water extraction distributions are shown in Table 4-3. The
percentage of water extracted from the top layer was greater for the
tree with grass cover.
Evapotranspiration data for the two grass cover treatments are
shown in Figures 4-12 and 4-13. Figure 4-12 shows ET rates of the 20
kPa grass covered lysimeter and Figure 4-13 shows ET rates of the 40 kPa
grass covered lysimeter. Comparison of Figure 4-12 and 4-13 shows that
the ET of the 20 kPa treatment was higher than that of the 40 kPa
treatment. The grass had more water available for ET in the 20 kPa
treatment.
Table 4-4 shows the water extraction distributions for the two
grass treatments. Both distributions were similar with both having
approximately the same soil water extraction percentages for each layer.
Model Verification
The accuracy of the numerical models in simulating the infiltration
and redistribution of soil water was determined by comparison with other
computer simulations from previous works in the literature.
The one-dimensional computer model was compared to the simulations
of Rubin and Steinhardt (1963) and Hiler and Bhuiyan (1971). Their work
provided data from soils with widly different hydraulic properties and
also provided simulation results which were used for model verification.
I
I
I
I
JULIAN DAY
Figure 4-10. Evapotranspiration rate from the 1986 field experiment with a
20 kPa soil water potential treatment young citrus tree with
grass cover.
Figure 4-11. Evapotranspiration rate from the 1986 field experiment with
a 20 kPa soil water potential treatment young citrus tree with
bare soil.
ro
73
Table 4-3. Two-dimensional distribution of water extraction for a 4
day drying cycle for young citrus trees with and without
grass cover.
Depth Depth Radius
Increment (mm)
(mm) (rrm)
200 600 200 600
Tree with Tree with
Grass Cover No Grass Cover
Relative Water Extraction
(%)
0-225
225
28.2
31.2
19.9
16.95
225-375
150
12.2
15.2
12.3
17.9
375-525
150
4.45
1.2
6.85
3.75
525-750
225
4.4
1.2
8.4
8.7
750-1050
300
0.75
1.2
2.6
2.65
Water Extraction
(nm)
0-225
225
1.6
1.8
0.9
0.8
225-375
150
0.7
0.9
0.6
0.8
375-525
150
0.3
0.07
0.3
0.2
525-750
225
0.3
0.07
0.4
0.4
750-1050
300
0.04
0.07
0.1
0.1
Water Extraction
(mm/mm)
0-225
225
0.0073
0.0081
0.0041
0.0035
225-375
150
0.0048
0.0059
0.0038
0.0057
375-525
150
0.0017
0.00046
0.0021
0.0011
525-750
225
0.0011
0.00031
0.0017
0.0018
750-1050
300
0.0001
0.00023
0.0004
0.0004
4
Figure 4-12.
JULIAN DAY
Evapotranspiration rate from the 1986 field experiment with
a 20 kPa soil water potential treatment with a grass cover.
-P»
4
3-|
Lü
170 171 172 173
JULIAN DAY
174
Figure 4-13. Evapotranspiration rate from the 1986 field experment with
a 40 kPa soil water potential treatment with a grass cover.
^1
CT»
76
Table 4-4. One-dimensional distribution of water extraction for a 4 dry
drying cycle for a grass cover crop at water depletion
levels of 20 kPa and 40 kPa.
Depth
Depth
Water
Water
Relative
Increment
Extraction
Extraction
Water
Extraction
(mm)
(mm)
(rim)
(mm/mm)
(*)
20 kPa Treatment Tree with
Grass Cover
0-225
225
6.3
0.028
57.9
225-375
150
2.0
0.014
18.6
375-525
150
0.9
0.006
8.3
525-750
225
1.1
0.005
9.9
750-1050
300
0.6
0.002
5.3
40 kPa Treatment Tree with
Grass Cover
0-225
225
3.9
0.0170
63.1
225-375
150
0.8
0.0055
13.5
375-525
150
0.5
0.0036
8.7
525-750
225
0.4
0.0019
6.9
750-1050
300
0.5
0.0016
7.8
77
Rubin and Steinhardt (1963) used an implicit solution of the
Richards equation in terms of water content. Their model was used to
study constant intensity rainfall infiltration on a Rehovot sand. The
soil data were presented as analytic functions which were tabulated for
use in this work. The soil hydraulic characteristics are shown in
Figures 4-14 and 4-15.
The simulated soil profile had a uniform initial soil water content
of 0.005. A uniform grid size of 10 mm was used. Figures 4-16 and 4-17
contain plots of soil water content profiles for infiltration into the
Rehovot sand. The figures show both the results from Rubin and
Steinhardt (1963) and this work. Figure 4-16 contains the soil water
content profiles for a constant infiltration rate of 12.7 mm/hr. Figure
4-17 contains the soil water content profiles for a constant
infiltration rate of 47 mm/hr. Results of this work are in excellent
agreement with those of Rubin and Steinhardt (1963).
The work of Hiler and Bhuiyan (1971) was also used to verify the
accuracy of the model developed in this work. They used a computer
model written in CMSP to solve the Richards equation. Surface
infiltration was simulated for two soils, Yolo light clay, and Adelanto
loam. The hydraulic characteristics of these soils are presented in
Figures 4-18 through 4-21.
The simulation results for these soils are presented in Figures 4-
22 through 4-24. Figures 4-22 and 4-23 show water content profiles with
time for the Yolo light clay at different initial conditions. Figure 4-
24 shows the soil water content profiles with time for the Adelanto
loam.
78
Comparison of results for the two soils showed excellent agreement
between Hiler and Bhuiyan's CMSP model and this work. The minor
differences between the results of the models could have resulted from
the interpolation of the soil properties.
The two-dimensional model was tested by simulating infiltration for
a Nahal Sinai sand (Bresler et al. 1971). Results were compared with a
steady-state analytical solution developed by Wooding (1968) and with an
ADI-Newton numerical simulation method presented by Brandt et al.
(1971).
These models simulated the infiltration of water into the soil from
a point source. The boundary condition used for the simulations was a
constant flux for the innermost surface grid (r = 0). The flux was set
equal to the saturated hydraulic conductivity of the soil. No flow
boundaries were used for the lower and radial boundaries. The simulated
soil profile had a uniform initial soil water content of 0.037. Figures
4-25 and 4-26 show the soil characteristics for the Nahal Sinai sand.
Figures 4-27 and 4-28 show a comparison between this work and that
of Wooding (1968) and Brandt et al. (1971). The data are presented in
dimensionless form with the relative water content defined as
S(0)/S(6sat) where S(0) was defined by Philip (1968) as
0
S(e) = \ D de (67)
0n
where D = the hydraulic diffusivity and 9n = the residual soil water
content. The total infiltration time was 57 min in Figure 4-27 and 297
min in Figure 4-28.
Figures 4-27 and 4-28 demonstrate that this model accurately
simulates the results of Brandt et al. (1971). Differences between the
SOIL WATER POTENTIAL (mm)
79
Figure 4-14. Soil water potential-soil water content relationship
for Rehovot sand.
HYDRAULIC CONDUCTIVITY (mm/hr)
80
Figure 4-15. Hydraulic conductivity-soil water content relationship
for Rehovot sand.
DEPTH (mm)
SOIL WATER CONTENT (mm/mm)
Figure 4-16. Simulated results of soil water content profiles for infiltration into a
Rehovot sand under constant rain intensity of 12.7 mm/hr.
DEPTH (mm)
Figure
SOIL WATER CONTENT (mm/mm)
17. Simulated results of soil water content profiles for infiltration into
Rehovot sand under constant rain intensity of 47 mm/hr.
SOIL WATER POTENTIAL (mm)
83
SOIL WATER CONTENT (mm/mm)
Figure 4-18. Soil water potential-soil water content relationship
for Yolo light clay.
HYDRAULIC CONDUCTIVITY (mm/hr)
Figure 4-19. Hydraulic conductivity-soil water content relationship
for Yolo light clay.
SOIL WATER POTENTIAL (mm)
85
Figure 4-20. Soil water potential-soil water content relationship
for Adelanto loam.
HYDRAULIC CONDUCTIVITY (mm/hr)
86
Figure 4-21. Hydraulic conductivity-soil water content relationship
for Adelanto loam.
DEPTH (mm)
0.0 0.1 0.2 0.3 0.4 0.5
SOIL WATER CONTENT (mm/mm)
Figure 4-22. Simulated results of soil water content profiles for infiltration into
a Yolo light clay with initial pressure potential at -66 kPa.
co
^4
Pages
Missing
or
Unavailable
RELATIVE EVAPOTRANSPIRATION RATE
TIME OF DAY
T I
25
Figure 4-32. Distribution of evapotranspiration through the daylight hours.
LO
SOIL HATER POTENTIAL (mm)
100
Figure 4-33. Soil water potential-soi1 water content relationship
for Arredondo fine sand.
HYDRAULIC CONDUCTIVITY (mm/hr)
.101
Figure 4-34. Hydraulic conductivity-soil water potential relationship
for Arredondo fine sand.
102
conductivity function produced results with the numerical model which
were in agreement with those from the observed field data.
One-dimensional Model
Time Steps. A sensitivity analysis was conducted to evaluate
optimum grid sizes and timesteps. The time step used was variable, and
the criteria for determining the variable time step were presented in
equation (31). Operation of the model provided data from which limits
for maximum and minimum time steps were determined. The maximum time
step for model operation was determined to be 1 hour. Operation of the
model with greater maximum time step produced varing results while
operation with smaller maximum time steps gave results in agreement with
the 1 hour timestep. The maximum time step was a upper limit for the
calculated time step. The maximum time step was used when water
movement in the soil profile was small. When water movements in the
soil profile were large, as during infiltration, a minimum time step was
determined. The minimum time step for the range of infiltration used in
the model was determined to be 0.0001 hours. The use of a smaller time
step slowed the model execution while not improving the mass balance.
Larger time steps caused the mass balance to increase, thus indicating a
reduction in model accuracy.
Grid Size. The optimum grid sizes were also determined by
sensitivity analysis. A grid size of 30 mm was determined to be
adequate to describe water movement in the soil profile for the one¬
dimensional model. Larger grid sizes decreased model execution time,
but the soil water distribution was not adequately described. Smaller
grid sizes increased the execution time while not improving the model
accuracy.
103
Model Operation. The accuracy of the one-dimensional model was
verified by simulating 35 days of observed water extraction,
infiltration and redistribution under a young citrus tree with grass
cover in a field lysimeter system. Initial conditions for the model
were obtained directly from the observed soil water potentials and
calculated root activity distributions. Observed average daily ET
values were input to the numerical model. A no flux boundary condition
was specified for the lower boundary. The surface boundary was set to a
constant flux during irrigation and no flux at the end of the
irrigation. Evapotranspiration from the top grid in the model was
accounted for in the soil water extraction term.
Simulations of the observed field data were in good agreement with
the experimental data. The simulated water contents were compared to
field data in Figures 4-35 to 4-39. These data represent the five soil
depths monitored with the data acquisition system. The simulated water
contents for the first drying cycle were in excellent agreement with the
field data.
An irrigation of 28 mm of water was applied to the lysimeter. The
model simulated this irrigation, redistribution and extraction for a
second drying cycle. The simulation results were in excellent agreement
with the field data for the second drying cycle also. The water
contents for the 450, 600 and 900 mm depths were in excellent agreement
with the observed field water contents throughout the drying cycle.
The water contents for the 150 mm depth was in good agreement throughout
most of the second drying cycle. The simulated water contents for the
300 mm depth were greater than the observed field data.
SOIL WATER CONTENT (mm/mm)
0.20
0.15
0.10
0.05
0.00
230 240 250 260 270
JULIAN DAY
Figure 4-35. Simulation results for the one-dimensional model with field data fo
the 150 mm depth.
— Field data
■ Model results
i • 1 • r
SOIL WATER CONTENT (mm/mm)
Figure 4-36. Simulation results for the one-dimensional model with field data for
the 300 mm depth.
SOIL WATER CONTENT (mm/mm)
0.2ft
— Field data
■ Model results
0.15-
0.10-
0.05.
230 240 250 260 270
JULIAN DAY
Figure 4-37. Simulation results for the one-dimensional model with field data
for the 450 mm depth.
O
O'»
SOIL WATER CONTENT (mm/mm)
Figure 4-38. Simulation results for the one-dimensional model with field data
for the 600 mm depth.
SOIL WATER CONTENT (mm/mm)
Figure 4-39. Simulation results for the one-dimensional model with field data
for the 900 mm depth.
109
The differences between the simulated and observed water contents
may have been due to the variations in the soil properties under field
conditions. The model assumed homogeneous and isotropic soil
conditions. This cause is supported by the random variations between
observed and simulated water contents. Overall the model simulated
changes in soil water movement and extraction with excellent agreement.
Two-dimensional Model
Time Steps. A sensitivity analysis was conducted to evaluate
optimum grid sizes and time steps for the two-dimensional model. The
time step used was variable, and the criteria for determining the
variable time step were presented in equation (31). Operation of the
model provided data from which limits for maximum and minimum time steps
were determined. These maximum and minimum time steps were determined
to be the same as used with the one-dimensional model. The maximum
allowable time step for model operation was determined to be 1 hour.
The minimum allowable time step for the range of infiltration used in
the model was determined to be 0.0001 hours.
Grid Size. The grid sizes for the two-dimensional model were also
determined by sensitivity analysis. A grid size of 30 mm was determined
to be adequate to accurately describe water movement in the vertical
direction. The radial grid size of 50 mm was determined to be adequate
to accurately describe water movement in the radial direction.
Model Operation. The two-dimensional model was used to simulate 35
days of water extraction, infiltration and redistribution under a young
citrus tree with grass cover in a field lysimeter. The initial
conditions for the model were determined directly from the observed soil
water potentials and root distributions. Observed daily ET amounts were
110
input to the model. Boundary conditions for the two-dimensional model
were all no flow boundaries when there was no infiltration. During
infiltration, a surface flux equal to the irrigation application rate
was used. When irrigation was completed, the surface boundary condition
was reset to a no flow condition.
Simulations of the observed field data were in excellent agreement
with the experimental data for the first drying cycle. The simulated
water contents were compared to the field data in Figures 4-40 to 4-44.
These data represent the 15 sites in the soil profile which were
monitored with the data acquisiton system. The observed and simulated
water contents for the drying cycle are presented for each of the four
radii measured. Simulated water contents were in excellent agreement
with the observed field water contents for almost all of the observed
sites during the first drying cycle. The 450 mm depth water contents
did not agree as well at the 100 mm radius for the later part of the
drying cycle.
Simulated soil water movement and extraction for the second drying
cycle were in good agreement with the observed field data. Simulated
water contents for the 100 and 300 mm radii were in excellent agreement
with field water contents for all five observed depths. The simulated
water contents for the outer radii were different than observed water
contents for several locations. These differences may be attributed to
variations in the soil properties or nonuniform application of
irrigation water at the surface. The soil properties were assumed to be
homogenous and isotropic. The model assumed that the irrigation
application was uniform. Overall the model results were in good
agreement with the observed field data.
SOIL WATER CONTENT (mm/mm)
0.20-
Field data
Model results
0.15-
0.15-
0.10-
0.05-
0.00.-
230
JULIAN DAY
Figure 4-40. Simulation results for the two-dimensional model with field data for
the 150 mm depth.
SOIL WATER CONTENT (mm/mm)
0.2a
0.15-
0.10-
0.05-
0.00- —
0.15-
0.10L
0.05-
o.oo4—
230
Figure 4-41.
Field data
Model results
JULIAN DAY
Simulation results for the two-dimensional model with field data for
the 300 mm depth.
SOIL WATER CONTENT (mm/mm)
JULIAN DAY
Figure 4-42. Simulation results for the two-dimensional model with field data for the
450 mm depth.
SOIL WATER CONTENT (mm/mm)
Figure 4-43. Simulation results for the two-dimensional model with field data for
the 600 mm depth.
SOIL WATER CONTENT (mm/mm)
Figure
0.20-
- - ■ — ' ■■—— —
Field data
Model results
—4,
¿ irr. i .■.■.■.■•4»
100 mm Radius
300 mm Radius
■ i i
^ U
i i l
500 mm Radius
700 mm Radius
230 240 250 260 22
< i |
tO 240 250 260 270
0.00-
0.00
JULIAN DAY
4-44. Simulation results for the two-dimensional model with field data for
the 900 mm depth.
116
Model Applications
The numerical models developed in this research can provide
researchers with numerical models of soil water movement to study
infiltration, water extraction and solute movement. The models may be
used to determine optimum irrigation scheduling for a wide variety of
soils and atmospheric conditions. The numerical models allow the
input of historic rainfall and ET data.
To demonstrate potential applications of the numerical models
developed, the one-dimensional numerical model was used to compare three
irrigation strategies. The first strategy studied was to schedule
irrigations to occur when the water potential in the soil profile
dropped below -10 kPa at a specified depth in the soil profile. The
second was to schedule irrigations to occur when the water potential in
the soil profile dropped below -20 kPa and the third when the water
potential dropped below -30 kPa. The numerical model was modified to
monitor specific depths in the soil profile to represent tensiometer (or
other soil water status monitoring equipment) locations as in field
operations. Depths of 300 and 900 mm were monitored in this study. An
irrigation depth of 25 mm was applied when the soil water potential at
any monitoring depth dropped below the set water potential. Rainfall
interactions were not considered in these simulations but could be taken
into account with the numerical model by inputting historical data.
The inputs to the model were the initial soil water potentials
which were assumed to be at field capacity. Evapotranspiration was
input on a daily interval, and for these simulations was assumed to be a
constant 4 mm/day. The surface boundary condition for the one¬
dimensional model was a no flux when irrigation was not occurring.
117
During irrigation the surface boundary condition was a flux boundary
with the flux equal to the irrigation application rate. The lower
boundary condition for the model was a gravity flow condition to
simulate deep percolation from the soil profile.
Simulation results for the three irrigation strategies were shown
in Figures 4-45 to 4-49. Figure 4-45 and 4-46 show the soil water
contents monitored at 300 and 900 mm. From the soil water contents, it
appears that the irrigations for the -10 kPa treatment were controlled
by the potentials at the 300 mm depth because of the frequent and small
irrigations applications. The irrigations for the -20 and -30 kPa
treatments appear to be controlled by the 900 mm depth because of the
less frequent and larger irrigation applications. These demonstrate
another use of the numerical model to determine the optimum depths to
monitor for different soils and different allowable soil water
depletions.
Figure 4-47 shows the soil water storage in the profile during the
simulation period. The soil water storage for the -10 kPa treatment was
consistently more uniform, while the soil water storage for the -20 and
-30 kPa treatments varied widely.
Figures 4-48 and 4-49 shows the irrigation applied and the drainage
out of the profile. The same amount of irrigation was applied for all
three treatments. Irrigation was applied more frequently for the -10
kPa treatment. Irrigations were less frequent for the -20 and -30 kPa
treatments but they required several irrigations to bring the soil water
potentials in the soil profile above the set limits. Figure 4-49 showed
the drainage loss out of the bottom of the profile for the different
SOIL WATER CONTENT (mm/mm)
0.20-
10 kPa
20 kPa
0 10 20 30 40 50 60
DAYS
Figure 4-45. Soil water contents for the three irrigation treatments at
the 300 mm depth.
00
SOIL WATER CONTENT (mm/mm)
Figure 4-46. Soil water contents for the three irrigation treatments for
the 900 mm depth.
‘119
STORAGE (mm)
Figure 4-47. Soil water storage for the three irrigation treatments.
ro
o
250
10 kPa
20 kPa
30 kPa
I 1 1 1 1 >
30 40 50 60
DAYS
Figure 4-48. Cumulative irrigation for the three irrigation treatments.
CUMULATIVE DRAINAGE (mm)
150
125
10 kPa
20 kPa
30 kPa
100.
75 -
/
Figure 4-49.
Cumulative drainage from the soil profile for the three
irrigation treatments.
123
treatments. The drainage out of the all of the treatments was large at
first because of the initial conditions. The drainage for the -10 kPa
teatment was small for the duration of the simulation period. Drainage
from the other treatments were higher because of the larger amounts of
water applied at each irrigation.
CHAPTER V
SUMMARY AND CONCLUSIONS
A data collection system was developed which provides
researchers with information needed for developing and validating models
of soil water movement, crop water use and evapotranspiration. Two
numerical models were developed to simulate soil water movement and
extraction. These numerical models may be used to compare irrigation
scheduling strategies.
Microcomputer-based Data Acquisition System
A low-cost microcomputer-based data acquisition system to
continuously record soil water potentials was developed and tested.
Excellent agreement was obtained between soil water potentials read with
the data acquisition system and those read manually using mercury
manometers. Typical agreement was within 0.47 kPa. The maximum,
deviation observed for all evaluations conducted was 1.76 kPa. The
microcomputer based data acquisition system has the ability to monitor
inputs and make decisions based upon the input data to control events in
the field.
Field Data
The data acquisition system developed was used to measure soil
water extraction patterns under young citrus trees and under a grass
cover crop grown in field lysimeters. Observations were monitored
during two drying cycles and one irrigation. The data collected were
evaluated and parameters determined for input to the numerical
simulation models. Evapotranspiration, soil water extraction rates and
124
125
relative soil water extractions were calculated from the soil water
potential measurements.
Modeling Soil Water Movement and Extraction
A one-dimensional numerical model was developed to simulate soil
water movement and extraction under a uniformly irrigated soil profile.
The model was verified using data from the literature. It was then used
to simulate the soil water movement and extraction from a young citrus
tree in a field lysimeter. Two drying cycles and one irrigation were
simulated. Simulated results were in excellent agreement with the data
collected using the microcomputer based data acquisition system.
A two-dimensional numerical model was developed to simulate soil
water movement and extraction under a nonuniformly irrigated soil
profile. The model was compared to both one- and two-dimensional models
of soil water movement from the literature. It was then used to
simulate the two-dimensional movement and extraction of soil water in a
field lysimeter with a young citrus tree. Two drying cycles and one
irrigation were simulated. The model results were in good agreement
with the data collected using the microcomputer based data acquisition
system.
Model Applications
The numerical models may be used to study irrigation scheduling
problems. A comparison of three different irrigation scheduling
strategies was conducted to demonstrate one use of the models. These
numerical models will also be useful in other applications such as the
modeling of solute movement in soils, modeling of fertilizer movement
under a trickle emitter, or as inputs to a larger crop modeling system.
APPENDIX A
Alternating Direction Implicit Finite Differencing
An alternating direction implicit solution method was developed to
solve the two-dimensional radial flow model. The solution method
divides the solution into two half cycles, thus two cycles are required
to solve for one time step. Equation (8) was solved in the radial
direction during the first half cycle and then solved in the vertical
direction during the second half cycle. The ADI solution for equation
(8) may be written as
k+1/2
k+1/2
k+1/2
k+1/2 8 h
1
1 3 k+1/2 3 h
1 3 k+1/2 3 h
C — =
- (
(r K _
) + (r K
—) )
3 t
r
2 3r 3 r
2 3 r
3 r
k+1
k
13 k+1 3 h
13 k 3 h
+
(
(K —)
+ (K _)
)
2 3 Z 3 z
2 3 Z 3 Z
k+1/2
k+1 k
1 3 K 1
3 K 3 K
+
(
+ — (
+ )
- S(z,r,t)
2 3 z 4
3 Z 3 Z
(68)
The finite difference equations used for the solution of equation (8)
were expressed as
k+1 k
(h - h )
i,j iJ
C
i,j
At
(69)
126
127
3 3h
—(K )
3 Z 3 Z
K ( h - h ) - K ( h - h )
Í+1/2J i+l,j i,j i-1/2,j i,j i-l,j
(70)
2
Az
13 3 h 1
_(_(rK _)) = (r K (h -h )
r 3r 3r 2 i,j+l/2 i,j+l/2 i,j+l i,j
r Ar
i.j
- r K ( h - h ) ) (71)
iJ-1/2 i,j-1/2 i.j i,j-1
K - K
3 K i-1/2, j 1+1/2,j
(72)
Equation (68) was divided into two parts and arranged such that the
unknowns are on the left side of the equations and the knowns are on the
right side such that
1 At
( 2
2 Ar
ri,j-l/2
r
i.j
k+1/2
K
i-1/2,j k+1/2
) h
i. J-1
C
i.j
k+1/2
1 At ( r K
i, j+1/2 iJ+1/2
+ (
k+1/2
+ r K )
i,j-1/2 1,j-1/2
+ 1 )
k+1/2
h
i.j
2
2
Ar r C
i.j i.j
128
k+1/2
1 At r K
i,j+1/2 i, j+1/2
2
2 a r r C
i,j i,j
k+1/2
h
k k
1 At K - K
k i-1/2, j i+1/2,j
= h +-— ( )+ (
i.j
2 A z C 4 A z
k+1/2 k+1/2
1 At K - K
i-1/2,j i+1/2,j
-)
1 A t K
+ ( -
i-1/2,j
2 A z C
k k
.) ( h - h )
1-1. j i.j
i,J
1 A t K
+ (
2 AZ C
i+1/2, j k k
) ( h - h )
1+1. J i.J
i.J
1 At
2 C
i.J
i.J
and
1 At K
( -
k+1
i-1/2, j k+1
) h
2 A z C
i.J
1-1, j
k+1 k+1
1 At ( K + K )
i-1/2,j i+1/2,j
+ (
2
AZ C
i,j
k+1
+ 1 ) h
i.j
(73)
2
129
1
- ( -
2
k
At K
i+1/2,j k+1
) h
2 1+1,j
Az C
i.j
1 At
+ _— (
4 Az
k+1/2 k+1/2
K - k
t-1/2, j i+1/2, j
)
C
i,j
1 A t
2 C
i.j
k
S
i.j
k+1/2
1
At r
i.j-1/2
K
I.j-1/2
2
2
A r r
C
i.j
i.j
1
At r
k+1/2
K
i.j+1/2
i.j+1/2
2
2 Ar r C
i.j 1
k+1/2 k+1/2
h - h )
i.j-1 i.j
k+1/2 k+1/2
h - h ) (74)
i.j+1 i.j
where i represents the vertical distance, j represents the radial
distance and k represents the time increment.
The system of equations which are formed by applying equations (73)
and (74) to each grid point produces a tridiagonal matrix at each half
time step. The tridiagonal matrices were solved using a Gauss
elimination method for tridiagonal matrices. Equations (73) and (74)
may be simplified such that equation (73) may be written as
k+1/2 k+1/2 k+1/2
- AR h + BR h - CR h = DR
i.j-1 i.j i.j+l i.j
(75)
130
where
k+1/2
1 A t r K
i.j-1/2 i.J-1/2
AR = ( )
2
2 a r r
i,J
i.J
(76)
k+1/2
1 At r K
l.j+1/2 i,j+1/2
CR = ( )
2
2 a r r
i.J
i.J
(77)
and
BR = 1 + AR + CR
(78)
1 At K - K
DR = h + (
i,j
2 AZ
i-1/2 ,j i+1/2,j
k+1/2 k+1/2
1 At K - K
i-1/2, j i+1/2,j
-) + (
4 AZ
1,J
-)
1 At K
+ (-
i-1/2,j
2 az C
k k
■) ( h - h )
1-1J l.j
i.J
1 At K
+ ( -
i+1/2, j
2 AZ C
i.J
1 At
k k
) ( h - h )
i+l»j 1.j
2 C
i.J
i.J
(79)
and equation (74) may be written as
131
- AV h
where
k+1
k+1
k+1
+ BV h
- CV h
1-1. j
I.j
1+1 .j
k+1
1
At K
AV =
(
1-1/2,j
2
2
A Z C
)
i,j
DV
i.J
(80)
(81)
k
1 At K
i+1/2J
C V = ( ) (82)
2
2 az C
i J
and
BV = 1 + AV + C V
(83)
DV
1 At
k
h + — (
i.j
4 AZ
k+1/2 k+1/2
K - K
1-1/2,j i+1/2, j
)
C ..
i,j
1 At
k
S
i.j
1 At r
i.j-1/2
2
2 A r r
i,j
k+1/2
K
i.j-1/2 k+1/2 k+1/2
) ( h - h )
i.j-1 1. J
C
I.j
k+1/2
1 At r K
i.j+1/2 i.j+1/2
+ ( )
2
2 Ar r C
i,j i.j
k+1/2 k+1/2
( h - h )
i,j+l 1,j
(84)
132
Boundary Conditions. The boundary condition for the surface layers
was a flux boundary condition which accounted for irrigation and
rainfall at the soil surface. At times when no irrigation occured the
surface boundary was a no flux boundary. Evaporation from the soil
surface was included in the water extraction term. The modification of
equations (75) and (80) to describe the surface boundary condition was
k+1/2 k+1/2 k+1/2
- AR h + BR h - CR h
l.j-1 l.j l.j+1
DR (85)
l.j
where
k+1/2
lit Qs - K 1 At Qs
k j 1+1/2,j j
DR = h + ( ) + (
l.j
2 az C 4 az
1J
- K
1+1/2,j
-)
l.j
1 At K
+ ( -
1+1/2,j k
2 az C
l.j
1 At
) ( h - hk ) Sk
1+1, j l.j l.j
2 C
l.j
(86)
and
k+1
BV h
l.j
where
k+1
CV h = DV
2,j l.j
BV = 1 + CV
(87)
(88)
and
133
DV =
k
h
l.j
1 At
+ __ (
4 Az
k+1/2
Qs - K
j 1+1/2,j
)
C
l.j
1 At
2 C
1J
k
S
l,j
k+1/2
1
At r
l.j-1/2
K
l.j-1/2
2
2
A r r
C
l,j
1, j
1
At r
k+1/2
K
1J+1/2
1, j+1/2
2
2
A r r
C
l.j
l.j
k+1/2 k+1/2
h - h )
l.j-1 l.j
k+1/2 k+1/2
h - h ) (89)
l.j+l l.j
The lower boundary condition was also represented as a flux
boundary. The lower boundary in this research was an impermeable
boundary which represented the bottom of the lysimeter. The
modifications of equations (75) and (80) to describe the lower boundary
condition were
k+1/2 k+1/2 k+1/2
- AR h + BR h - CR h = DR
n,j-l n,j n,j+1 n,j
where AR, BR and DR are as previously described and
DR
1 At
k
h +
n,j
2 az
k
K
n-1/2J
Qb 1 At
j
) + - — (
k+1/2
K - Qb
n-1/2, j j
)
C 4 Az C
n,j n,j
(90)
134
k
1 At
n-l/2,j k k k
) ( h - h ) S
2 n-1, j n,j n,j
2 az C 2 C
n»j n,j
1 At K
+ (
(91)
and
k+1 k+1
- AV h + BV h
n-1, j n,j
where AV was previously described and
BV = 1 + AV
and
k+1/2
1 At K - Qb
k n-1/2,j j
DV = h + (
n.j
4 az C
n,j
DV
n,J
1 At
-) S
2 C
n»J
(92)
(93)
n,J
k+1/2
1 At r K
n,j-1/2 n,j-1/2
2
2 at r
k+1/2 k+1/2
) ( h - h )
n,j-l n,j
n,J
n,J
k+1/2
1 At r K
n.j+1/2 n,j+1/2
2
2 a r r C
k+1/2 k+1/2
) ( h - h )
n,j+l n, j
n,J
n,J
(94)
The boundary conditions for the radial boundaries were both no flux
boundaries. The no flux boundary condition for the outer radius (j=m)
135
represented the outer wall of the lysimeter and equations (75) and (80)
were written as
k+1/2 k+1/2
- AR h + BR h = DR (95)
i,m-l i,m i,m
where
and
A V h
BR = 1 +
AR
(96)
k+1
k+1
k+1
+
BV h
- CV h = DV
(97)
i-1 ,m
i,m
i+l,m i,m
where
k+1/2 k+1/2
1 At K - K 1 At
k i-l/2,m i+l/2,m k
DV = h + ( ) S
i,m i ,m
4 A z C 2 C
i,m i,m
k+1/2
1 At r K
i.m-1/2 i,m-l/2 k+1/2 k+1/2
+ ( ) ( h - h )
2 i,m-l i,m
2 a r r
i ,m
i ,m
(98)
The boundary condition for the inner radius was also represented as
a no flux boundary. For the inner radius r(i,j) = 0, therefore its
inclusion into equation (8) would result in a division by zero.
Therefore, the second term in equation (8) was rewritten as described in
equations (50) and (51) such that the finite difference equations for
the ADI solution of equation (51) were expressed as
136
where
and
where
DV =
k+1/2 k+1/2
BR h - CR h = DR
i ,1 i ,2 i ,1
k+1/2
1 At K
i+1/2,1 k+1/2
CR = ( ) h
2 1,2
2 Ar C
1,1
BR = 1 + CR
k+1 k+1 k+1
AV h + BV h - CV h = DV
1-1,1 i ,1 1+1,1 1,1
k+1/2 k+1/2
1 At K - K
k i -1/2,1 i+1/2,1
h + - — ( )
M
4 az C
1,1
1 At
2 C
1,1
1,1
(99)
(100)
(101)
(102)
k+1/2
1
At K
i,1+1/2
k+1/2 k+1/2
+ (-
)
( h - h
2
1,2 1,1
2
Ar C
i,l
(103)
Updating of the soil parameters, surface infiltration, mass balance
and time steps were implimented as described for the one-dimensional
soil water flow model
APPENDIX B
Listing of Soil Water Movement and Extraction Model
137
138
C
c
C This program models the water movement and extraction in both a one-
C and a two-dimensional soil profiles. The model was designed to
C study water movement and extraction and how they affect irrigation
C scheduling. The program used an input file to input the parameters
C for the model. The input file is interactively asked for during
C program execution.
C
C The soils used in this model are assumed to homogenous and isotropic
C The surface boundry is a flux boundary which may be specified in the
C input file. The boundary is assumed to be a no flux unless an
C irrigation is occuring. The lower boundary may be either a no flux
C boundary or a gravity flow. The radial boundaries for the two-
C dimensional model are both no flow boundaries.
C
C
C CONT
C COND
C DTIRR
C DZ
C DEPTH
C DR
C DT
C DTMAX
C DTMIN
C ENDTIME
C HEAD
C HEDT
C ILOW
C I CONST
C ISOIL
C MBFILE
C NL
C NR
C R
C QS
C S
C SWC
C SWCT
C TIME
C TOUT
C TOD
C WC
C WCNT
C
PARAMETER (IZ=50,1R=16,1B=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
INTEGER I, NUM, I OUT,ISTOP,I COUNT
CHARACTER*64 MBFILE
COMMON /PARA/WC(IZ,IR,3),C0ND(IZ,IR,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),I CONST,NR,DR(IR) ,R(IR),IL0W
COMMON/TIM/DT,TIME .TOUT(90) .ENDTIME ,DTMAX ,DTMIN,TOD,ITIME ,QM,DTIRR
COMMON /TABLES/WCNT(120),C0NT(120),SWCT(120) ,HEDT(120)
= table of conductivities
= the hydraulic conductivity
= the duration of the irrigation event
= the vertical grid increment
= the depth of the grids
= the radial grid increment
= the time step
= the maximum time step
= the minimum time step
= the ending simulation time
= the matric water potential
= table of soil matric potentials
= the status of the lower boundary condition
= the status of the surface boundary condition
= the soil type
= the output masbalance file
= the number of vertical grids
= the number of radial increments
= the radius of the grids
= the surface flux
= the soil water extraction
= the specific water capacity
= table of specific water capacities
= the current simulation time
= the times as which the outputs should occur
= the time of day
= The water content
= table of water contents
o o o
139
$ ,NDATA,MBFILE, ISO II-
CO MMON /MASB2/DIFF,PSTOR.FLOWIN,EXTRACT,STORAGE.FLOWOUT
COMMON /WCC/WKEEP,FLOWl(IZ,IR)
COMMON /CONV/DIF,ITER,MAXITER
COMMON /QS1/RDF(IZ,IR),WCWP,WCFC,ICF
COMMON /ET1/TSNR(90),TSNS(90),ETP(90),RAIN(90),DURJSTOP
COMMON /DUMY/IFIEL0(25) ,JFIELD(25) ,IFS,RTIME,1 OPT
$ ,ITEN,ITF(25),JTF(25),THRES,RRATE,ISCH(90)
CALL INPUT
CALL MASBAL(O)
WRITE(*,*)' THE INITIAL MASS BALANCE IS '
WRITE (*,90) TIME , DT,DIFF,PSTOR .FLOWOUT
WRITE(11,*)' THE INITIAL MASS BALANCE IS '
WRITE (11,90) TIME ,DT,DIFF,PSTOR, FLOWOUT
WRITE(*,91)STORAGE,FLOWIN,EXTRACT
WRITE(11,91)STORAGE,FLOWIN,EX TRAC T
TNEW = RTIME
IF(IOPT.EQ.0)WRITE(15,924)TNEW,(WC(IFIELD(I),JFIELD(I),1),1=1,IFS)
START OF SIMULATION LOOP
IT IME=1
MITER = 0
IDAY = TIME / 24.0 + 1.0
88 CONTINUE
CALL FLUX(ETPFLX, ETP(IDAY),TSNR(IDAY),TSNS(IDAY),TOD)
CALL QSS( ETPFLX)
CALL IRRIGATE( IDAY, IRCODE )
C ONE-DIMENSIONAL MODEL
IF(NR.EQ.l) THEN
CALL TSTEP1D( IRCODE )
CALL REDST1D
END IF
C TWO-DINENSIONAL MODEL
IF(NR.GE.3) THEN
CALL TSTEP2D( IRCODE )
CALL REDST2D
END IF
C
TIME = TINE + DT
IDAY = TIME / 24.0 +1.0
TOD = TOD + DT
IF(T0D.GE.24.0) TOD = TOD - 24.0
C
IF(TIME.EQ.TOUT(ITIME).OR.TIME.EQ.ENDTIME) THEN
C OUTPUT LOOP
ITIME =ITIME+ 1
WRITE(*,90)TINE ,DT,DIFF,PSTOR,FLOWOUT
WRITE(11,90)TIME,DT,DIFF,PSTOR,FLOWOUT
90 FO RMAT (/,
$' TINE1 ,T25, '= 1 , F12.5 ,' (HRS) DT',9X,'= \F12.5
(HRS)',/,
$’ MASS BALANCE --> ERROR = \E12.5,' (mm3) TERROR
$ E12.5,/,' DRAINAGE = ' ,E12.5)
140
C
WRITE(*,91)STORAGE,FLOWIN.EXTRACT
WRITE(11,91)STORAGE.FLOWIN.EXTRACT
91 FORMAT(1 STORAGE = \E12.5,' CUM_INFIL = '.E12.5,
$' CUM_EXTRACT = ' ,E12.5,/)
C
WRITE(11,9000)TINE
9000 FORMATC TIME = ' ,F12.5)
DO 12 1=1, NL
WRITE(*,11) DEPTH(I) ,(WC(I,J,1) ,J=1,NR)
IF(IOPT.EQ.l)WRITE(15,11)(WC(I,J ,1),J=NR,1,-l),(WC(I,J ,1) ,J=1,NR)
12 WRITE( 11,11) DEPTH(I) ,(WC(I,J,1) ,J=1,NR)
11 F0RMAT(10E12.5)
95 FORMAT(E10.4 ,10X,E10.4,10X,E10.4)
C
TNEW = RTINE + TINE / 24.0
IF(IOPT.EQ.0)WRITE(15,924)TNEW,(WC(IFIELD(I),JFIELD(I),1),I=1,IFS)
IF(I0PT.EQ.3)WRITE(15,924)TNEW,(WC(IFIELD(I),JFIELD(I),1),I=1,IFS)
$ .STORAGE, FLOWIN, EXTRACT, FLOWOÜT
924 FORMAT(18X,F9.4,20E12.5)
END IF
C
c
MITER = MITER + 1
IF(MITER/5*5.EQ.MITER) THEN
OPEN(12,FILE='TIME.0UT1',IOSTAT=IC)
IF(IC.NE.O) GO TO 9008
WRITE(12,9005)TIME,MI TER
9005 FORMATC TIME = ' ,F12.5,' ITER = ',110)
CLOSE(12)
9008 END IF
C
IF(ENDTIME.LE.TINE) GO TO 99
GO TO 88
99 CONTINUE
OP EN(17,FILE='FINAL.POT')
DO 66 1=1 ,NL
66 WRI TE (17 ,11) ( HEAD(I,J,1),J=1,NR)
CLQSE(17)
CLOSE(ll)
END
C -
C SUBROUTINE TO READ IN THE INPUT VALUES FOR THE SIMULATION
C AND INITIALIZE THE VARIABLES
C
SUBROUTINE INPUT
PARAMETER (IZ=50,IR=16,IB=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
CHARACTER*64 INPFILE,MBFILE,OUTFILE,WDFILE,DATAFILE
CHARACTER*! Z
REAL IWCONT.ISTOR, HE(6,5),EX(6,5)
INTEGER 11(6),JJ(5)
COMMON /TABLES/WCNT(120),C0NT(120),SWCT(120),HEDT(120)
$ ,NDATA,MBFILE,ISO IL
141
COMMON /PARA/WC(IZ,IR,3),COND(IZ,IR ,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL,DZ(IZ),DEPTH( IZ),QS( IR) ,S(IZ,IR),ICONST,NR,DR(IR),R(IR),ILOW
COMMOM/irM/DT.TIME,T0UT(90).ENDTIME,DTMAX,DTMIN,T0D,ITIME,QM,DTIRR
COMMON /MASB/ IWCONT( IZ, IR) ,AREA( IR) .TAREA,ISTOR
COMMON /DUMY/IFIELD(25) ,JFIELD(25) ,IFS,RTIME,IOPT
$ ,ITEN,ITF(25),JTF(25),THRES,RRATE,ISCH(90)
COMMON /CONV/DIF,ITER,MAX ITER
COWON /QS1/RDF(IZ,IR) ,WCWP,WCFC,ICF
COMMON /ET1/TSNR(90),TSNS(90),ETP(90),RAIN(90),DUR,TSTOP
COMMON /WCC/WKEEP,FL0W1(IZ,IR)
PIE = 3.1415926
WRITE(*,*)1 ENTER THE INPUT DATA FILE1
READ(*,1Q00) DATAFILE
0PEN(8,FILE=DATAFILE)
READ(8 ,901) Z
READ(8,1005) ISOIL.INPFILE
1000 F0RMAT(A64)
1005 FORMAT( II ,1X ,A64)
OPEN( 9, FI LE= INP FILE)
C
C READ IN AND WRITE THE SOILS DATA
C
READ(9,11) NDATA
11 FORMAT(15)
READ(9,09)(HEDT(I),WCNT(I),SWCT(I),CONT(I),1=1,NDATA)
09 F0RMAT(4E12.5)
900 FORMAT(' ')
CLOSE!9)
c
READ(8,901) Z
READ(8,1000)MBFILE
READ!8,1000)WDFILE
0PEN(11,FILE=MBFILE)
READ!8,901) Z
READ(8,' (II ,1X ,A64)') IOPT.OUTFILE
0PEN(15,FILE=OUTFILE)
READ(8,*) IFS
READ( 8,*) RTIME , (I FI ELD( I) ,JFIELD( I) , 1=1,IFS)
READ(8,901) Z
READ(8,*)TII€ .ENDTIME ,TOD
READ(8,901)Z
READ(8,*)WKEEP,WCWP,WCFC
READ(8,901) Z
READ(8,*)TOUT(1)J0UT(2),T0UT(3),T0UT(4)
IF(TOUT(2).NE.-99.0) GO TO 9000
DO 7000 I = 1,90
IF(T0UT(1)*FL0AT(I).GT.ENDTIME) GO TO 7000
TOUT(I) = TOUT(1)*FLOAT(I)
7000 CONTINUE
9000 CONTINUE
READ(8,901) Z
READ(8,*)DTMAX,DTMIN,MAXITER,DTIRR,DUR
WRITE(*,710)MBFILE.TIME,TOUT(1),ENOTIME,DTMAX,DTMIN,MAX ITER
WRITE(11,710)MBFILE,TIME,TOUT(1).ENDTIME,DTMAX,DTMIN,MAX ITER
142
710 FORMAT(/,* THE Í4ASS BALANCE DATA IS IN FILE ',A40,/,
$' THE STARTING SIMULATION TIME IS ' ,E12.5,/,
$' THE OUTPUTS WILL BE AT INTERVALS 1 ,E12 .5,/,
$' THE ENDING TIME WILL BE ',E12.5,/,
$' THE MAXIMUM TIME STEP IS ' ,E12.5,/,
$' THE MINIMUM TINE STEP IS ' ,E12.5,/,
$' THE NUMBER OF ITERATIONS ',112,/)
C
C READ IN THE NUMBER OF GRIDS AND THE GRID SIZES
C
READ(8 ,901) Z
READ(8,*) NL ,DZ(1), NR, DR(1)
WRITE(*,720)NL,DZ(1),NR,DR(1)
720 FORMAT(1 THE NUMBER OF VERTICAL INCREMENTS ',110,/,
$ ' THE VERICAL SPACING IS ’,E12.5,/,
$ ' THE NUMBER OF RADIAL INCREMENTS ',110,/,
$ ' THE RADIAL SPACING IS ' ,E12.5,/)
C
DEPTH(l) = - DZ(1) / 2
DO 200 I = 2,NL
DZ(I) = DZ(I-l)
DEPTH(I) = DEPTH(1-1) - ( DZ(I) + DZ(I-l) ) / 2
200 CONTINUE
C
IF(NR.EQ.2) NR = 1
AREA(l) = 1.0
TAREA = 1.0
IF(NR.GE.3) THEN
R(1) = DR(1) / 2
AREA(l) = PIE * (R(l)+DR(l)/2)**2
DO 210 I = 2, NR
DR(I) = DR(I-l)
R(I) = R(I-l) + ( DR(I) + DR(I-l) ) / 2
AREA( I)= PIE* ( (R(I)+DR(I)/2)**2 -(R( I)-DR( I)/2)**2)
210 CONTINUE
TAREA = PIE * ( DR(1) * NR )**2
WRITE (* ,735)
735 FORMAT(/,1 THE RADIUS AND AREAS ',/)
WR ITE (* ,740) (I ,R( I) ,AREA( I) ,1=1 ,NR)
740 FORMATC RADIUS(1 ,13,') = 1 .E12.5 ,10X,' AREA = ‘ ,E12.5)
END IF
C
C INPUT THE INITIAL POTENTIALS FOR THE SYSTEM
C
READ(8,901)Z
C READ(8,910) (( HEAD(I,J,1),J=1,IR),I=1,IZ)
C
READ(8,*) NLN.NRN.INV
READ(8,*) ( 11( I), 1=1 ,NLN )
READ(8,*) ( JJ(J),J=1,NRN )
DO 1010 1=1,NLN
1010 READ(8,*) ( HE(I,J),J=1,NRN)
READ(8 ,901)Z
C
C READ IN THE INFILTRATION RATE
C
C READ(8,910) (QS(J),J=1,IR)
READ(8,901) Z
READ(8,901)Z
READ(8,*) ICONST
READ(8,901)Z
C
C READ IN THE CONDITION OF THE LOWER BOUNDRY CONDITION
C AND THE STATUS OF THE CORRECTION FACTOR
C
READ(8,*)ILOW, ICF
READ(8,901) Z
C READ IN THE ROOT DISTRIBUTION FUNCTION
C
C RE AD (8,910) ((RDF( I ,J) ,J=1 ,IR) , 1=1 ,IZ)
DO 1011 1=1,NLN
1011 READ(8,*) ( EX(I,J),J=1,NRN)
C
DO 811 K=1,3
DO 811 1=1,IZ
DO 811 J=1,1R
HEAD(I ,J ,K)=0.0
RDF(I,J) = 0.0
S(I,J) = 0.0
811 CONTINUE
C
N = 0
DO 400 1=1,NLN
DO 400 L=1,11(I)
N = N + 1
M = 0
DO 400 J=1 ,NRN
DO 400 K=1 ,JJ(J)
M = M + 1
HEAD(N,M,1) = HE(I ,J)
HEAD(N,M,2) = HE(I,J)
HEAD(N,M,3) = HE(I,J)
K1 = M
IF(INV.EQ.l) K1 = NR + 1 - M
AAA = FL0AT(K1**2 - (Kl-1)**2) / FLOAT( NR**2)
IF(NR.EQ.l) AAA = 1.0
RDF(N,M) = EX(I ,J) / 11( I) * AAA
400 CONTINUE
C
910 F0RMAT(5F10.0)
901 FORMAT(Al)
775 FORMAT( 10E12.5)
C
CALL UPDATE (1)
ISTOR =0.0
DO 310 1=1,NL
DO 310 J=1 ,NR
IWCONT(I ,J) =WC( I ,J ,1)
144
ISTOR = ISTOR + WC(I,J ,1)*DZ(I)*AREA(J)
310 CONTINUE
C
110 FORMAT( 17 ,3X,4(El2.5 ,3X) )
READ(8,*) ITEN, THRES, RRATE, ISC
RE AD (8,*) ( ITF(I) ,JTF( I) ,1=1, ITEN )
DO 1015 1=1,IDAY
1015 ISCH(I) = ISC
CL0SE(8)
C
C READ IN WEATHER DATA
C
0PEN(10, FILE=WDFILE)
IDAY = ENDTIME / 24.0 + 1.0
REA0( 10 ,1001 ,END=1002) ( ETP(I) ,TSNR(I) ,TSNS(I) ,RAIN(I) ,1=1,IDAY)
1002 WRITE(11,1003)
WRITE(11,1001)( ETP(I),TSNR(I),TSNS(I),RAIN(I),1=1,IDAY)
1001 F0RMAT(4F10.2)
CLOSE(IO)
1003 FORMAT(/,1 THE WEATHER DATA ',/)
C
RETURN
END
C
SUBROUTINE FLUX (ETPFLX , ETP, TSNR ,TSNS ,TOD)
C
C PROGRAM TO CALCULATE THE DISTRIBUTION OF POTENTIAL
C ET DURING THE DAY. THE DISTRIBUTION IS ASSUMED TO
C BE A SINE FUNCTION.
C -
C DEFINITIONS
C TSNR = TIME OF SUNRISE
C TSNS = TIME OF SUNSET
C ETP = POTENTIAL ET MM
C ETPRMX = MAXIMUM ET RATE MM/HR
C DAYL = DAY LENGTH HOURS
C ETPFLX = CURRENT ET FLUX MM/HR
C
IMPLICIT REAL (A-H.O-Z)
PI = 3.141592654
DAYL = TSNS - TSNR
ETPRMX = 2 * ETP / DAYL
OMEGA = 2 * PI / DAYL
TD = TOD - TSNR
TI = TOD
IF(TI .LT.TSNR.OR.TI.GT.TSNS) TD = 0.0
ETPFLX = (ETPRMX/2)*SIN((0MEGA*TD) - PI/2.) + ETPRMX/2
END
C
C SUBROUTINE TO CALCULATE THE WATER UPTAKE RATE FOR THE
C ONEDIMENSIONAL SOIL PROFILE
C
SUBROUTINE QSS(TRANS)
PAR/4 METER (IZ=50,IR=16,IB=IR*2+1 ,MAX=IZ*IR)
145
IMPLICIT REAL (A-H.O-Z)
REAL TRANS,RASW(IZ,IR),A(IZ.IR),K1,ISTOR,IWCONT
COMMON /PARA/WC(IZ,IR,3) ,COND(IZ,IR,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R(IR),I LOW
COMMON /MASB/IWCONT(IZ,IR),AREA(IR),TAREA,ISTOR
COMMON /QS1/RDF(IZ,IR),WCWP,WCFC,ICF
C
IF(TRANS.EQ.O.O) THEN
DO 30 1= 1,NL
DO 30 J= 1,NR
30 S(I,J) = 0.0
RETURN
END IF
K1 = 10.0
SUM = 0.0
DO 10 1 = 1,NL
DO 10 J=1 ,NR
RASW( I ,J) = ( WC(I,J,1) - WCWP ) / ( WCFC - WCWP )
IF(RASW(I,J).GT.l.0) RASW(I ,J) = 1.0
IF(RASW(I,J).LE.O.0) THEN
RASW(I,J) =0.0
A(I,J) = 0.0
GO TO 100
END IF
A(I, J) = TRANS * TAREA * RASW(I,J) ** ( TRANS /(RASW(I ,J)*K1))
100 S(I,J) = A( I ,J) * RDF (I ,J)
SUM = SUM + S(I,J)
10 CONTINUE
C
IF( SUM.EQ .0.0 .OR. ICF.EQ .1) THEN
CF = 1.0
ELSE
CF = TRANS * TAREA / SUM
END IF
SUM2 = 0.0
DO 20 1=1 ,NL
DO 20 J=1 ,NR
S(I,J) = S(I,J) / DZ(I) / AREA(J) * CF
SUM2 = SUM2 + S(I,J)*DZ(I)*AREA(J)
20 CONTINUE
RETURN
END
C
SUBROUTINE TRID(NG,IF,A,B,C,D,V)
C
C THIS SUBROUTINE SOLVES THE TRIDIAGONAL COEFFICIENT MATRIX
C WHICH RESULTS FROM A SET OF SIMULTANEOUS EQUATIONS WHICH
C CAN BE PUT INTO THE FOLLOWING FORM
C - A*V1 + B*V2 - C*V3 = D
C
PARAMETER (IZ=50)
IMPLICIT REAL (A-H.O-Z)
DIMENSION A(IZ),B(IZ),C(IZ),D(IZ),V(IZ),BETA(IZ),GAMMA(IZ)
IFP1 = IF+1
146
C NGM1 = NG-1
C -
C COMPUTE INTERMEDIATE ARRAYS BETA AND GAMMA
C
BETA(IF) = B(IF)
GAMMA(IF) = D(IF)/BETA(IF)
DO 100 I = IFPl.NG
BETA(I) = B(I)-A(I)*C(I-1)/BETA(I-1)
GAMMA(I) = (D(I)+A(I)*GAMMA(I-1))/BETA(I)
100 CONTINUE
C
C COMPUTE FINAL SOLUTION VECTOR V
C
V(NG) = GAMMA (NG)
LAST = NG-IF
DO 200 K = 1 ,LAST
I = NG-K
V(I) = GAMMA(I)+C(I)*V(1+1)/BETA(I)
200 CONTINUE
RETURN
END
C
C SUBROUTINE TO CALCULATE A MASS BALANCE OF THE SYSTEM SIMULATED
C
SUBROUTINE MASBAL(IRESET)
PARAMETER (IZ=50,IR=16,IB=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
REAL IWCONT,ISTORAGE
COMMON /PARA/WC(IZ,IR,3),COND(IZ,IR,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),I CONST,NR,DR(IR),R(IR),I LOW
COMMON/TIM/DT,TIME,T0UT(90).ENDTIME,DTMAX,DTMIN,TOD,ITIME ,QM,DTIRR
COMMON /MASB/ IWCONT(IZ.IR),AREA(IR)JAREA,ISTORAGE
COMMON /MASB2/DIFF,PST0RAGE,FL0WIN,EXTRACT,STORAGE,FLOWOUT
C
C INITIALIZE THE CUMULATIVE INFILTRATION AND EXTRACTION
C
IF(IRESET.EQ.0) THEN
FLOWIN = 0.0
FLOWOUT =0.0
EXTRACT =0.0
DIFF=0 .0
STORAGES.0
DO 20 J=1 ,NR
DO 20 1=1,NL
20 STORAGE = STORAGE + WC(I ,J ,1)*DZ(I)*AREA(J)
DIFF = STORAGE - ISTORAGE
PSTORAGE = DIFF / ISTORAGE * 100.0
RETURN
END IF
DIFF= 0.0
STORAGES.0
C
DO 10 J=1,NR
FLOWIN = FLOWIN + DT*QS(J)*AREA(J)
147
IF(ILOW.EQ.1) FLOWOUT = FLOWOUT + DT*(COND(NL,J,1) )
DO 10 1=1,NL
STORAGE = STORAGE + WC(I,J,1)*DZ(I)*AREA(J)
EXTRACT = EXTRACT + S(I,J)*DZ(I)*AREA(J)*DT
10 CONTINUE
DIFF = STORAGE - (ISTORAGE + FLOWIN - EXTRACT -FLOWOUT )
PSTORAGE = DIFF / (ISTORAGE + FLOWIN - EXTRACT - FLOWOUT) * 100.0
RETURN
END
C
SUBROUTINE IRRIGATE( IDAY, IRCODE )
C
C SUBROUTINE TO CALCULATE THE AMOUNT OF IRRIGATION REQUIRED
C FOR THE MODEL BASED ON THE MAXIMUM SOIL WATER POTENTIAL ALLOWED
C
PARAMETER (IZ=50,IR=16,IB=IR*2+1 ,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
COMMON /PARA/WC(IZ,IR,3),COND(IZ,IR,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ.IR),ICONST,NR,DR(IR),R(IR),1 LOW
COMMON/TIM/DTJIME ,T0UT(90) ,ENDTIME .DTMAX ,DTMIN ,TOD ,ITIME ,QM,DTIRR
COMMON /ET1/TSNR(90) ,TSNS(90) ,ETP(90) ,RAIN(90) ,DUR,TSTOP
COMMON /QS1/RDF(IZ,IR),WCWP,WCFC,ICF
COMMON /DUMY/1FIELD(25),JFIELD(25),IFS ,RTIME,IOPT
$ ,ITEN,ITF(25),JTF(25),THRES,RRATE,ISCH(90)
COMMON /MASB/ IWCONT(IZ,IR),AREA(IR),TAREA,ISTORAGE
C
AMOUNT =0.0
DEFICIT = 0.0
IF(ISCH(IDAY).EQ.O) THEN
DO 10 J=1 ,NR
DO 10 1=1,NL
DEFICIT = DEFICIT + ( WCFC - WC(I,J,1) ) * DZ(I) * AREA(J)
10 CONTINUE
AMOUNT = DEFICIT / TAREA
DO 20 1= 1 ,ITEN
IF (ABS(HEAD(ITF(I),JTF(I),1)).GE.ABS(THRES)) THEN
RAIN(IDAY)=AMOUNT
ISCH(IDAY)= 1
END IF
20 CONTINUE
END IF
C
IF( RAIN(IDAY).LE.0.0) THEN
DO 100 1=1,NR
QS( I) = 0.0
100 CONTINUE
IRCODE = 0
RETURN
END IF
C
IF( IRCODE.EQ.O) THEN
TSTART = TIME
C DUR = RAIN(IDAY) / RRATE
RAIN(IDAY) = RRATE
non nonnn
148
TSTOP = TSTART + DUR
C
IRCOOE = 1
00 200 1=1 ,NR
QS( I) = RAIN( IDAY )
200 CONTINUE
END IF
C
IF( TINE .GE. TSTOP ) THEN
IRC0DE = 0
DO 300 1=1,NR
QS( I) = 0.0
300 CONTINUE
RAIN(IDAY) =0.0
END IF
RETURN
END
SUBROUTINE TO CALCULATE THE TIME STEP BASED ON THE
MAXIMUM ALLOWABLE CHANGE IN WATER CONTENT
FOR THE ONE-DIMENSIONAL MODEL
SUBROUTINE TSTEP1D( IRCODE )
PARAMETER (IZ=50 ,IR=16,IB=IR*2+1 ,MAX=IR*IZ)
REAL K1,K2
COMMON/PARA/WC(IZ,IR,3) ,C0ND( IZ,IR,3),SWC(IZ,IR,3) ,HEAD( IZ,IR,3),
$NL,DZ(IZ) ,DEPTH(IZ),QS(IR),S(IZ,IR),I CONST,NR,DR(IR),R(IR),1 LOW
COMMON/TIM/DT.TIME,T0UT(90),ENDTINE ,DTMAX,OTMIN ,TOD,ITIME ,QM,DTIRR
COMMON /WCC/WKEEP, FL0W1( IZ ,IR)
COMMON /ET1/TSNR(90),TSNS(90),ETP(90),RAIN(90),DUR,TSTOP
CALCULATE THE MAXIMUM FLUX FROM THE INFILTRATION AND EXTRACTION
RATE = ABS(QS(1))
EXT = ABS(S(1,1)*DZ(1))
DO 10 1=1,NL
IF(EXT.LT.ABS( S( 1,1 )*DZ( I))) EXT = ABS( S( 1,1) *DZ( I) )
10 CONTINUE
C
IF(RATE.LT.EXT) RATE = EXT
QM = RATE
C
DO 100 1=1,NL
FLOW1 (1,1)=0.0
C
IF(I.EQ.l) THEN
FLOWl(I.l) = QS( 1)
$ + ( COND( 1,1,1) + COND( 1+1,1,1) ) / 2
$ *((HEAD(1+1,1,1)-HEAD(1,1,1))/DZ(I)-1.0)
GO TO 900
END IF
C
IF(I.EQ.NL) THEN
FLOWl(I.l) =
149
$ -1.0*( CONO(1,1.1) + COND(1-1,1,1) ) / 2
$ *((HEAD(1,1 ,1)-HEAD(I-1,1,1))/DZ(I)-1.0)
IF (I LOW. EQ. 1) FLOWl(I.l) = FLOWl(I.l) - COND( 1,1,1)
GO TO 900
END IF
C
FL0W1 (1,1) =
$ -1.0*( COND( 1,1,1) + COND( 1-1,1,1) ) / 2
$ *((HEAD(1,1,1)-HEADC1-1,1,1) )/DZ( I)-1.0)
$ + ( COND( 1,1,1) + COND( 1+1,1,1) ) / 2
$ *( (HEAD( 1+1,1,1) -HEAD( 1,1,1) )/DZ( I)-1.0)
C
900 IF(ABS(FLOW1 (1,1)) .GT.QM) GO TO 200
GO TO 100
200 QM = ABS(FL0W1 (1,1))
100 CONTINUE
C
IF(QM.LE.O.O) GO TO 400
DT = WKEEP * DZ(1) / QM
GO TO 500
400 DT = DTMAX
500 IF(DT.GT.DTMAX) DT = DTMAX
IF(DT.LT.DTMIN) DT = DTMIN
IF(IRCODE.EQ.1) DT = AMIN1( DTIRR, DT )
IF(DT+TIME.GT.T0UT(ITIME).AND.TOUT(ITIME).GT.O .0)
$ DT = TOUT(ITIME) - TINE
IF(DT+TINE.GT.ENDTIME) DT = ENDTINE - TINE
IF(DT+TINE.GT.TSTOP.AND.IRCODE.EQ.1) DT = TSTOP - TIME
IF(DT.LT.O.O) DT = DTMIN
C WRITE(*,191) DT,QM
191 FORMAT(/,****** DT = *,F10.8 , ’ QM = ',E12.5)
RETURN
END
C
C SUBROUTINE TO CALCULATE THE NEW WATER POTENTIALS
C FOR THE ONE-DINENSIONAL MODEL
C
SUBROUTINE REDST1D
PARAMETER (IZ=50,IR=16 ,IB=IR*2+1 ,MAX=IR*IZ)
IMPLICIT REAL (A-H.O-Z)
REAL K1,K2,D(IZ),A(IZ),B(IZ),C(IZ),V(IZ)
REAL CONTAB(120),WCCTAB(120) ,SPCTAB(120),PRTAB(120)
CHARACTER*64 MBFILE
COMMON /TABLES/WCCTAB,CONTAB,SPCTAB,PRTAB,NDATA,MBFILE,ISOIL
COMMON/PARA/WC(IZ,IR,3),COND( IZ,IR,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R(IR),ILOW
COMMON/TIM/DT,TIME,T0UT(90),ENDTIME,DTMAX,DTMIN,TOD,ITIME,QM,DTIRR
CONMON /CONV/DIF,I TER,MAX ITER
C
C SET UP THE INTERIOR GRID POINTS
C
DO 10 KL=2 ,NL-1
Kl=( COND(KL,1,1) + C0ND(KL-1,1,1) ) / 2
K2=( COND( KL+1,1,1) + COND( KL ,1,1) ) / 2
150
Cl = SWC(KL,1,1)
C0NST1 = DT / DZ(KL)**2/Cl
C
A(KL) = C0NST1 * K1
C(KL> = C0NST1 * K2
B(KL) = 1.0 + A(KL) + C(KL)
D(KL) = HEAD(KL,1,1) + DT/DZ(KL)/C1*( K1 - K2 )
$ - DT * S(KL,1)/Cl
10 CONTINUE
C
C SET UP THE BOUNDRY CONDITIONS
C -
IF(ICOPiST.EQ.2) GO TO 800
C
C NO FLUX BOUNDRY CONDITION AT THE SURFACE
C
KL = 1
K2=( CO ND(KL+1,1,1) + COND(KL ,1,1) ) / 2
K1 = 0.0
Cl = SWC (KL ,1,1)
C0NST1 = DT / DZ(KL)**2/C1
C
A(KL) = 0.0
C(KL) = C0NST1 * K2
B( KL) = 1.0+ C (KL)
D(KL) = HEAD(KL,1,1) + DT/DZ(KL)/Cl*( QS(1) - K2 )
$ - DT * S(KL ,1)/C1
GO TO 900
C
C CONSTANT POTENTIAL BOUNDRY CONDITION AT THE SURFACE
C -
800 KL = 1
K2=( CO'ND( 2,1,1) + COND( 1,1,1) ) / 2
K1 = C0ND( 1,1,1)
Cl = SWC(1,1,1)
C0NST1 = DT / DZ(1)**2/C1
C
A(l) = C0NST1 * K1
C(l) = C0NST1 * K2
B (1) = 1.0 + C(l)
D( 1) = HEAD(1,1,1) + DT/DZ(1)/C1*( K1 - K2 )
$ - DT * S(1,1)/C1
C
C NO FLUX BOUNDRY CONDITION AT LOWER BOUNDRY
C -
900 KL = NL
Kl=( COND(KL ,1,1) + COND( KL-1,1,1) ) / 2
K2= 0.0
IF( ILOW.EQ .1) K2 = COND( KL ,1 ,1)
Cl = SWC(KL ,1,1)
C0NST1 = DT / DZ(KL)**2/Cl
C
A(KL) = C0NST1 * K1
B(KL) = 1.0 + A(KL)
151
D(KL) = HEAD(KL,1,1) + DT/DZ(KL)/Cl*( K1 - K2 )
$ - DT * S(KL,1)/C1
C
C CALCULATE THE NEW VALUES FOR THE HEADS
C -
IF = 1
C IF(ICONST.EQ.2) IF = 2
CALL TRID(NL,IF,A,B,C,D ,V)
DO 700 KL = ICONST.NL
HEAD(KL ,1,3)=V(KL)
700 CONTINUE
C
IF(HEAD(1,1,3).GT.PRTAB(l).AND.ICONST.EQ .1) ICONST = 2
IF (ICONST.EQ.l) GO TO 2000
Q =0.0
Q = AMAX1( QS(1), QM )
C
2000 CONTINUE
CALL UPDATE (3)
DO 5000 1=1,NL
HEAD (1,1,1) =HEAD (1,1,3)
HEAD( 1,1,2)=HEAD( 1,1,3)
WC(I,1,1) =WC( 1,1,3)
WC(1,1,2) =WC(1,1,3)
COND( 1,1,1)=CQND( 1,1,3)
COND(1,1,2) = C0ND(1,1,3)
SWC( 1,1,1) = SWC (1,1,3)
SWC(1,1,2) =SWC(1,1,3)
5000 CONTINUE
CALL MASBAL(l)
C
RETURN
END
C
C SUBROUTINE TO CALCULATE THE TIME STEP BASED ON THE MAXIMUM
C NET FLOW RATE INTO EACH GRID
C FOR THE TWO-DIMENSIONAL MODEL
C
SUBROUTINE TSTEP2D( IRCODE )
PARAMETER (IZ=50,IR=16,IB=IR*2+1 ,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
REAL K1,K2,K4,K5,NFR
COMMON /WCC/WKEEP, FL0W1( IZ ,IR)
COMMON /ET1/TSNR(90),TSNS(90),ETP(90),RAIN(90),DUR,TSTOP
COMMON /PARA/WC(IZ,IR,3),C0ND(IZ,IR,3),SWC(IZ.IR,3),HEAD(IZ,IR,3),
$NL ,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R(IR),ILOW
COMMON/TIM/DT,TIME,T0UT(90),ENDTIME,DTMAX,DTMIN,TOD,ITIME,QM,DTIRR
NFR =0.0
DO 100 1=1,NL
DO 100 J =1,NR
IF(I.EQ.l) THEN
Q1 = QS(J)
ELSE
C
152
Q1 = -1.0 * (( COND(I,J,1) + C0ND(1-1,J ,1) ) / 2)
$ *((HEAD(I,J,1)-HEAD(1-1,J,1))/DZ(I)-1.0)
END IF
IF(I.EQ.NL) THEN
Q2 = 0.0
IF(ILOW.EQ.l) Q2 = COND( I,J,1)
ELSE
Q2= -1.0 * (( COND( 1+1 ,J ,1) + COND( I ,J ,1) ) / 2)
$ *( (HEAD( 1+1 ,J ,1) -HEAD( I ,J ,1) )/DZ( I)-1.0)
END IF
IF(J.EQ.l) THEN
Q3 = 0.0
ELSE
Q3= -1.0 * (( COND(I ,J ,1) + COND( I ,J-1,1) ) / 2)
$ *((HEAD(I ,J ,1) - HEAD(I ,J-1,1))/(R( J) -R( J-l)) )
END IF
IF(J .EQ.NR) THEN
Q4 = 0.0
ELSE
Q4 = -1.0 * (( CO ND( I, J ,1) + COND( I ,J+1,1))/ 2)
$ *( (HEAD( I ,J+1,1)-HEAD( I ,J ,1))/ (R( J+1) -R( J)) )
END IF
FL0W1 (I ,J) = (Q1 - Q2)/DZ( I) + (Q3 - Q4)/DR(J) - S( I ,J)
IF(ABS(FLOW1 (I ,J)) .GT.NFR) GO TO 200
GO TO 100
200 NFR = ABS(FL0W1(I ,J))
100 CONTINUE
C
IF(NFR.LE.O.O) GO TO 400
DT = WKEEP / NFR
GO TO 500
400 DT = DTMAX
500 IF(DT.GT.DTMAX) DT = DTMAX
IF(DT.LT.DTMIN) DT = DTMIN
IF(IRCODE.EQ.1) DT = AMIN1( DTIRR, DT )
IF( DT+TIME.GE.TOUT(ITIME).AND.TOUT(ITIME).GT.O.0)
$ DT= TOUT(ITIME) - TIME
IF(DT+TIME.GT.TSTOP.AND.IRCODE.EQ.1) DT = TSTOP - TINE
IF( DT+TIME.GE.ENDTIME) DT = ENDTIME - TIME
IF( DT.LE.O.0 ) DT = DTMIN
RETURN
END
C***
C
C ADI SOLUTION METHOD OF RICHARDS EQUATION
C
SUBROUTINE REDST2D
PARAMETER (IZ=50,1R=16,1B=IR*2+1 ,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
REAL K1,K2,K3,K4,V(IZ),A(IZ),B(IZ),C(IZ),D(IZ),HSTAR(IZ ,IR)
COMMON/PARA/WC(IZ,IR,3),COND(IZ ,IR,3),SWC(IZ,IR,3),HEAD(IZ,IR,3),
$NL ,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R( IR),I LOW
COMMON/TIM/DT,TIME,T0UT(90),ENDTIME,DTMAX,DTMIN,TOD,ITIME,QM,DTIRR
COMMON /CONV/DIF,I TER,MAX I TER
153
C
C SET UP THE INTERIOR GRID POINTS FOR THE FIRST HALF CYCLE
C SOLVE IN THE RADIAL DIRECTION
C
ITER=0
700 DO 20 KL = 1,NL
IF(KL.EQ.l) GO TO 25
IF(KL.EQ.ML)GO TO 35
DO 10 J=2,NR-1
I = J
Kl=( COND(KL,J,1) + C0ND(KL-1,J,1) ) / 2
' K2=( C0ND(KL+1 ,J ,1) + C0ND(KL,J,1) ) / 2
K4=( COND(KL ,J+1,2) + C0ND(KL,J,2) ) / 2
K3=( C0ND(KL,J,2) + COND(KL,J-l ,2) ) / 2
Cl=( SWC(KL,J,1) + SWC(KL ,J ,3) ) / 2
Rl=( R(J+l) + R(J) ) / 2
R2=( R(J) + R(J-l) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl / R(J)
A(I) = -1.0* C0NST2 * R2 * K3 / 2
C(I) = -1.0* C0NST2 * R1 * K4 / 2
B( I) = 1.0 + A(I) + C(I)
D( I) = HEAD(KL ,J ,1) + C0NST1 * (K2 + Kl) * HEAD(KL,J ,1) / 2
$ - C0NST1 * ( Kl * HEAD(KL-1,J ,1) + K2 * HEAD(KL+1 ,J ,1) ) / 2
$ + DT/DZ(KL)/C1*( Kl - K2 )/2 + DT/DZ(KL)/Cl*( Kl - K2 ) / 4
$ - DT/C1 * S(KL ,J) / 2
10 CONTINUE
C
C NO FLUX BOUNDRY CONDITION AT INNER RADIUS
C
J = 1
I = J
Kl = (COND(KL,J ,1)+C0ND(KL-1 ,J ,1))/ 2
K2=( C0ND(KL+1 ,J ,1) + C0ND(KL,J,1) ) / 2
K4=( COND(KL ,J+1,2) + C0ND(KL,J,2) ) / 2
K3= K4
Cl =( SWC(KL.J.l) + SWC(KL ,J ,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl
A( I) = 0.0
C(I) = -1.0* C0NST2 * K4 / 2
B( I) = 1.0 + C(I)
D( I) = HEAD(KL ,J ,1) + C0NST1*(K2+K1) * HEAD(KL,J ,1) / 2
$ - C0NSTl*Kl*HEAD(KL-l,J,l)/2 - C0NST1*K2*HEAD(KL+1,J ,l)/2
$ + DT/DZ(KL)/C1*( Kl - K2 )/2 + DT/DZ(KL)/Cl*( Kl - K2 ) / 4
$ - DT/Cl * S(KL,J) / 2
C
C
C NO FLUX BOUNDRY CONDITION AT OUTER RADIUS
C
J=NR
I = NR
Kl=( COND(KL ,J,1) + COND(KL-1 ,J ,1) ) / 2
K2= (CON D( KL+1 ,J ,1) + C0ND(KL,J,1) ) ' 2
154
K3=( COND(KL ,J ,2) + CQND(KL ,J-1,2) ) / 2
K4 = 0.0
Cl =( SWC(KL ,J ,1) + SWC(KL ,J ,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl
C
A(I) = -1.0* C0NST2 * K3 / 2
C(I) = 0.0
B( I) = 1.0 + A( I)
D( I) = HEAD(KL ,J ,1) + C0NST1*(K2+K1) * HEAD(KL ,J ,1) / 2
$ - C0NST1*K1*HEAD(KL-1,J,1)/2 - C0NST1*K2*HEAD(KL+1,J ,l)/2
$ + DT/DZ(KL)/Cl*( K1 - K2 )/2 + DT/DZ(KL)/Cl*( K1 - K2 )/4
$ - DT/C1 * S(KL ,J) / 2
GO TO 50
C
C NO FLUX BOUNDRY CONDITION AT THE SURFACE
C
25 00 30 J=1,NR
I = J
K2=( C0ND(KL+1 ,J ,1) + C0ND(KL,J,1) ) / 2
K1 = QS(J)
IF(J.NE.NR) K4=( CO ND (KL ,J+1,2) + C0ND(KL,J,2) ) / 2
IF(J.EQ.NR) K4 = 0.0
IF(J.NE.NR) R1 = ( R(J+1) + R(J) ) / 2
IF(J.NE.l) K3 = ( C0ND(KL ,J ,2) + C0ND(KL ,J-1,2) ) / 2
IF(J.EQ.l) K3 = 0.0
IF(J.NE.l) R2 = ( R(J) + R(J-l) ) / 2
Cl=( SWC(KL,J,1) + SWC(KL,J,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
IF(J.NE.l.AND.J.NE.NR) C0NST2 = -1.0 * DT / DR(J)**2 / Cl / R(J)
IF(J.EQ.l.OR.J.EQ.NR) C0NST2 = -1.0 * DT / DR(J)**2 / Cl
A(I) = -1.0* C0NST2 * R2 * K3 / 2
IF(J.EQ.l) A(I) = -1.0* C0NST2 * K3 / 2
C(I) = -1.0* C0NST2 * R1 * K4 / 2
IF(J.EQ.NR) C(I) = -1.0* C0NST2 * K4 / 2
B(I) = 1.0 + A(I) + C(I)
D(I) = HEAD(KL,J,1) + C0NST1 * K2 * HEAD(KL ,J ,1) / 2
$ - C0NST1 * K2 * HEAD(KL+1,J,1) / 2
$ + DT/DZ(KL)/C1*(K1 - K2)/2 + DT/DZ(KL)/Cl*(K1 - K2)/4
$ - DT/Cl * S(KL ,J) / 2
30 CONTINUE
GO TO 50
C
C NO FLUX BOUNDRY CONDITION AT LOWER BOUNDRY
C
35 DO 40 J=1 ,NR
I = J
Kl=( C0ND(KL,J,1) + C0ND(KL-1 ,J ,1) ) / 2
K2= 0.0
IF(J.NE.NR) K4=( C0ND(KL ,J+1,2) + C0ND(KL,J,2) ) / 2
IF(J.EQ.NR) K4=0 .0
IF(J.NE.NR) R1 = (R(J+l) + R(J) ) / 2
IF(J.NE.l) K3=( COND(KL ,J ,2) + COND(KL ,J-1,2) ) / 2
IF(J.EQ.l) K3 = 0.0
155
IF(J.NE.l) R2=( R(J) + R(J-l) ) / 2
Cl = ( SWC(KL,J,1) + SWC(KL,J,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
IF(J.NE.l.AND.J.NE.NR) C0NST2 = -1.0 * DT / DR(J)**2 / Cl / R(J)
IF(J.EQ.l.OR.J.EQ.NR) C0NST2 = -1.0 * DT / DR(J)**2 / Cl
A(I) = -1.0* C0NST2 * R2 * K3 / 2
IF(J.EQ.l) A(I) = -1.0* C0NST2 * K3 / 2
C(I) = -1.0* C0NST2 * R1 * K4 / 2
IF(J.EQ.NR) C(I) = -1.0* C0NST2 * K4 / 2
B( I) = 1.0 + A (I) + C( I)
D(I) = HEAD(KL ,J,1) + C0NST1 * K1 * HEAD(KL ,J ,1) / 2
$ - C0NST1 * K1 * HEAD(KL-1 ,J ,1) / 2
$ + DT/DZ(KL)/C1*(K1 - K2)/2 + DT/DZ(KL)/Cl*(K1 - K2)/4
$ - DT/C1 * S(KL,J) / 2
40 CONTINUE
****************************************
50 IF = 1
CALL TRID( NR,IF,A, B, C, D, V )
DO 55 J=1 ,NR
HEAD(KL ,J ,2) = V( J)
55 CONTINUE
20 CONTINUE
C
C SET UP THE INTERIOR GRID POINTS FOR THE SECOND HALF
C IN THE VERTICAL DIRECTION
C
DO 60 J=1 ,NR
IF(J.EQ.l) GO TO 12
IF(J.EQ.NR)GO TO 14
DO 70 KL=2,NL-1
I = KL
Kl=( COND(KL ,J ,3) + C0ND(KL-1,J ,3) ) / 2
K2=( C0ND(KL+1 ,J ,3) + C0ND(KL,J,3) ) / 2
K4=( C0ND(KL,J+1,2) + C0ND(KL,J,2) ) / 2
K3=( COND(KL,J ,2) + COND(KL,J-1,2) ) / 2
Rl=( R(J+l) + R(J) ) / 2
R2=( R(J) + R(J-l) ) / 2
Cl = ( SWC(KL ,J ,1) + SWC(KL ,J ,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl / R(J)
A(I) = -1.0* C0NST1 * K1 / 2
C(I) = -1.0* C0NST1 * K2 / 2
B(I) = 1.0 + A(I) + C(I)
D(I) = DT/DZ(KL)/Cl*( K1 - K2 ) / 4.0
$ - C0NST2*R2*K3*HEAD(KL,J-1,2)/2 - C0NST2*R1*K4*HEAD(KL,J+1,2)/2
$ + HEAD(KL,J,2) + C0NST2 * ( R2*K3 + R1*K4 ) * HEAD(KL ,J ,2) / 2
$ - DT/Cl * S(KL,J) / 2
70 CONTINUE
C
C SET UP THE BOUNDRY CONDITIONS
C
C NO FLUX BOUNDRY CONDITION AT THE SURFACE
C
KL = 1
o o o o o o
156
I = KL
K2=( C0ND(KL+1 ,J,3) + C0ND(KL,J,3) ) / 2
K1 = QS(J)
K4=( COND( KL ,J+1,2) + C0ND(KL,J,2) ) / 2
Rl=( R(J+l) + R(J) ) / 2
K3=( COND(KL ,J ,2) + COND(KL ,J-1,2) ) / 2
R2=( R(J) + R(J-l) ) / 2
C1 = ( SWC(KL,J,1) + SWC(KL ,J ,3) ) / 2
CONST1 = -1.0 * OT / DZ(KL)**2 / C1
C0NST2 = -1.0 * DT / DR(J)**2 / C1 / R(J)
A( I) = 0.0
C(I) = -1.0* C0NST1 * K2 / 2
B(I) = 1.0 + C(I)
D(I) = DT/DZ(KL)/C1*( K1 - K2 ) / 4
$ - C0NST2*R2*K3*HEAD(KL,J-1,2)/2 - C0NST2*R1*K4*HEAD(KL,J+1,2)/2
$ + HEAD(KL,J,2) + C0NST2 * ( R2*K3 + R1*K4 ) * HEAD(KL ,J,2) / 2
$ - DT/C1 * S(KL,J) / 2
NO FLUX BOUNDRY CONDITION AT LOWER BOUNDRY
KL = NL
I = KL
Kl=( COND(KL,J,3) + COND(KL-1 ,J ,3) ) / 2
K2= 0.0
K4=( COND(KL ,J+1,2) + COND( KL ,J ,2) ) / 2
R1 = (R(J+l) + R(J) ) / 2
K3=( COND(KL ,J ,2) + COND(KL ,J-1,2) ) / 2
R2=( R(J) + R(J-l) ) / 2
Cl = ( SWC(KL,J,1) + SWC(KL ,J,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl / R(J)
A( I) = -1.0* C0NST1 * K1 / 2
C( I) = 0.0
B(I) = 1.0 + A(I)
D(I) = DT/DZ(KL)/Cl*( K1 - K2 ) / 4
$ - C0NST2*R2*K3*HEAD(KL,J-1,2)/2 - C0NST2*R1*K4*HEAD(KL,J+l ,2)/2
$ + HEAD(KL ,J ,2) + C0NST2 * ( R2*K3 + R1*K4 ) * HEAD(KL,J ,2) / 2
$ - DT/C1 * S(KL ,J) / 2
GO TO 65
NO FLUX BOUNDRY CONDITION AT INNER RADIUS
12 DO 41 KL=1 ,NL
I = KL
IF(KL.NE.l) K1 = (COND(KL,J,3)+C0ND(KL-l ,J ,3))/ 2
IF(KL.EQ.l) K1 = QS(J)
IF(KL.NE . NL) K2= ( C0ND(KL+1 ,J ,3) + C0ND(KL,J,3) ) / 2
IF(KL.EQ.NL)K2 = 0.0
K4=( COND(KL,J+1,2) + C0ND(KL,J,2) ) / 2
K3= K4
Cl = ( SWC(KL,J,1) + SWC(KL ,J ,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl
A(I) = -1.0* C0NST1 * K1 / 2
157
IF(KL.EQ.l) A(I) = -1.0* 0.0
C(I) = -1.0* C0NST1 * K2 / 2
IF(KL.EQ.NL) C(I) = -1.0* 0.0
B (I) = 1.0 + A (I) + C( I)
D(I) = DT/DZ(KL)/Cl*( K1 - K2 ) / 4 - C0NST2*K4*HEAD(KL,J+1,2)/2
$ + HEAD(KL,J,2) + C0NST2 * K4 * HEAD(KL,J ,2) / 2
$ - DT/C1 * S(KL,J) / 2
41 CONTINUE
GO TO 65
C
C NO FLUX BOUNDRY CONDITION AT OUTER RADIUS
C
14 DO 61 KL=1,NL
I = KL
IF(KL.NE.l) Kl=( COND(KL ,J ,3) + COND(KL-1 ,J ,3) } / 2
IF(KL.EQ.l) Kl= QS(J)
IF(KL.NE.NL) K2=(C0ND(KL+1 ,J ,3) + C0ND(KL,J,3) ) / 2
IF(KL.EQ.NL)K2 = 0.0
K3=( COND(KL ,J ,2) + COND( KL ,J-1,2) ) / 2
K4 = 0.0
Cl = ( SWC(KL ,J ,1) + SWC(KL ,J ,3) ) / 2
C0NST1 = -1.0 * DT / DZ(KL)**2 / Cl
C0NST2 = -1.0 * DT / DR(J)**2 / Cl
A(I) = -1.0* C0NST1 * K1 / 2
IF(KL.EQ.l) A(I) = -1.0* 0.0
C(I) = -1.0* C0NST1 * K2 / 2
IF(KL.EQ.NL) C(I) = -1.0* 0.0
B(I) = 1.0 + A(I) + C(I)
DC I) = DT/DZ(KL)/Cl*(K1 - K2) /4 - C0NST2* K3* HEAD(KL,J-1,2)/2
$ + HEAD(KL ,J ,2) + C0NST2 * K3 * HEAD(KL ,J ,2) / 2
$ - DT/C1 * S(KL,J) / 2
61 CONTINUE
C
65 IF = 1
CALL TRID(NL,IF, A, B, C, D, V)
DO 71 KL=1,NL
HSTAR(KL,J)= V(KL)
71 CONTINUE
60 CONTINUE
C -
DIF =0.0
DO 506 1=1,NL
DO 506 J=1,NR
DIF = DIF + ABS( (HSTAR( I ,J) -HEAD( I ,J ,3) )/HEAD( I ,J ,3) )
HEAD( I ,J ,3)=HSTAR( I ,J)
506 CONTINUE
C
IF(DIF.LT.l .0E-6) GO TO 701
CALL UPDATE (2)
ITER= ITER + 1
IF(ITER.LT.MAX ITER) GO TO 7 00
C
701 CONTINUE
CALL UPDATE(3)
OC-JOOOOOO o o o
158
DO 80 1=1,NL
DO 80 J=1 ,NR
HEAD( I ,J ,1)=HEAD( I ,J ,3)
HEAD( I ,J ,2}=HEAD( I ,J ,3)
WC(I,J,1) =WC(I ,J ,3)
WC(I ,J ,2) =WC(I,J ,3)
C0ND( I ,J ,1)=C0ND( I ,J ,3)
C0ND(I ,J ,2) = C0ND(I ,J ,3)
SWC(I,J,1) =SWC( I ,J ,3)
SWC(I ,J ,2) =SWC(I ,J ,3)
80 CONTINUE
CALL MASBAL(l)
RETURN
END
ROUTINE TO UPDATE THE SOIL PROPERTIES
SUBROUTINE UPDATE ( IL00P )
PARAMETER (IZ=50,IR=16,IB=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
CHARACTER*64 MBFILE
REAL WCNT(120),C0NT(120) ,SWCT(120) ,HEDT( 120)
COMMON /TA BLES/WCNT,C0NT,SWCT,HEDT,NDIM,MBFILE,I SOIL
IF(ISOIL.EQ.l) CALL TBLARR ( ILOOP )
IF(ISOIL.EQ .2) CALL TBLREH ( ILOOP )
IF(ISOIL.EQ.3) CALL TBLOAM ( ILOOP )
IF(ISOIL.EQ .4) CALL TBLOAM ( ILOOP )
IF(ISO IL.EQ.5) CALL TBLYOL ( ILOOP )
IF(ISOIL.EQ .6) CALL LINEAR ( ILOOP )
RETURN
END
SUBROUTINE TBLARR ( ILOOP )
THIS SUBROUTINE USES AN INDEX ROUTING TECHNIQUE WITH
TABULATED DATA TO FIND THE CORRESPONDING VALUES OF
WATER CONTENT, SPECIFIC WATER CAPACITY AND HYDRAULIC
CONDUCTIVITY FOR A GIVEN VALUE OF SOIL MATRIC POTENTIAL
FOR ARREDONDO FINE SAND.
PARAMETER (IZ=50,IR=16,IB=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H,0-Z)
CHARACTER*64 MBFILE
REAL WCNT(120),CONT(120),SWCT(120),HEDT(120)
COMMON/PARA/WCN(IZ,IR,3),CON(IZ,IR,3),SWC(IZ,IR,3),HED(IZ,IR,3),
$ NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R(IR),ILOW
COMMON/TABLES/WCNT,CONT ,SWCT ,HEDT ,NDIM,MBFILE ,ISOIL
HEDSML = HEDT(118)
HEDLGE = HEDT(l)
DO 100 J=I LOOP ,3
DO 100 K=1 ,NR
DO 100 I = 1 ,NL
IF (HED (I,K ,J) .LE. HEDSML) HED(I,K,J) = HEDSML
o o o o o o o
159
IF (HED(I,K,J) .GE. HEDLGE) HED(I,K,J) = HEDLGE
RINDEX = AL 0G10(-1.*HED(I,K ,J))*4.E01-8.1EOl
INDEX = RINDEX
IF(INDEX.LT.l) INDEX = 1
IF (HED(I,K,J) .LT. HEDT(INDEX+1)) INDEX = INDEX+1
IF (HED(I,K,J) .GT. HEDT(INDEX)) INDEX = INDEX-1
IF(INDEX.LT.l) INDEX = 1
IF (INDEX .GE. 118) INDEX = 118
RM = RINDEX-INDEX
RATIO = AL0G10(HED(I,K,J)/HEDT(INDEX))/ALOG1O(HEDT(INDEX+1)/
+ HEDT(INDEX))
WCN( I ,K ,J) = WCNT( INDEX)+RATIO* (WCNT( INDEX+1)-WCNT( INDEX))
SWC(I ,K,J) = SWCT(INDEX)*1.E01**(RATI0*AL0G10(SWCT(INDEX+1)/
+ SWCT(INDEX)))
CON(I ,K,J) = CONT(INDEX)*1.E01**(RATI0*AL0G10(C0NT(INDEX+1)/
+ CONT(INDEX)))
100 CONTINUE
RETURN
END
C
SUBROUTINE TBLOAM ( ILOOP )
THIS SUBROUTINE USES AN INOEX ROUTING TECHNIQUE WITH
TABULATED DATA TO FIND THE CORRESPONDING VALUES OF
WATER CONTENT, SPECIFIC WATER CAPACITY AND HYDRAULIC
CONDUCTIVITY FOR A GIVEN VALUE OF SOIL MATRIC POTENTIAL
FOR EITHER OF THE LOAMS.
PARAMETER (IZ=50 ,IR=16 ,1 B=IR*2+1 ,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
CHARACTER*64 MBFILE
REAL WCNT(120),C0NT(120),SWCT(120),HEDT(120)
COMMON/PARA/WCN(IZ,IR,3),CON(IZ,IR,3),SWC(IZ,IR,3),HED(IZ,IR,3),
$ NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R(IR),ILOW
COMMON/TA BLES/WCNT,CONT,SWCT,HEDT.NDIM.MBFILE,ISOIL
HEDSML = HEDT(49)
HEDLGE = HEDT(l)
DO 100 J = I LOOP ,3
DO 100 K = 1,NR
DO 100 I = 1,NL
IF (HED(I,K,J) .LE. HEDSML) HED(I,K,J) = HEDSML
IF (HED( I ,K ,J) .GE. HEDLGE) HED(I,K,J) = HEDLGE
RINDEX = AL0G10(-1 .*HED( I ,K,J) )*1 .E01-2.0E01
INDEX = RINDEX
IF (HED(I,K,J) .LT. HEDT(INDEX+1)) INDEX = INDEX+1
IF (HED(I,K,J) .GT. HEDT(INDEX)) INDEX = INDEX-1
IF (INDEX .GE. 49) INDEX = 49
IF (INDEX .LE. 1) INDEX = 1
RM = RINDEX-INDEX
RATIO = AL0G10(HED( I, K,J)/HEDT( INDEX) )/AL0G10( HEDT( INDEX+1)/
+ HEDT(INDEX))
WCN(I.K.J) = WCNT( INDEX) +RATI0* (WCNT( INDEX+1)-WCNT( INDEX) )
SWC(I ,K ,J) = SWCT(INDEX)*!.EOl**(RATI0*AL0G10(SWCT(INDEX+1)/
160
+ SWCT(INDEX)))
CON(I,K,J) = C0NT(INDEX)*1,E01**(RAT10*AL 0G10(CO NT(INDEX+1)/
+ CONT{INDEX)})
100 CONTINUE
RETURN
END
C
SUBROUTINE TBLREH ( ILOOP )
C THIS SUBROUTINE USES AN INDEX ROUTING TECHNIQUE WITH
C TABULATED DATA TO FIND THE CORRESPONDING VALUES OF
C WATER CONTENT, SPECIFIC WATER CAPACITY AND HYDRAULIC
C CONDUCTIVITY FOR A GIVEN VALUE OF SOIL MATRIC POTENTIAL
C FOR REHOVOT SAND
C
PARAMETER (IZ=50,IR=16,IB=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
CHARACTERS MBFILE
REAL WCNT(120),CONT{120),SWCT(120),HEDT(120)
COMMON/PARA/WCN(IZ,IR,3),C0N(IZ,IR,3),SWC(IZ,IR,3),HED(IZ ,IR,3),
$ NL ,DZ( IZ) ,DEPTH( IZ) ,QS( IR) ,S(IZ,IR) , ICONST,NR,DR( IR) ,R(IR) ,ILOW
COMMON/TABLES/WCNT, CONT,SWCT,HEDT,NDIM,MBFILE,ISOIL
C
HEDSML = HEDT(59)
HEDLGE = HEDT(l)
DO 100 J = ILOOP ,3
DO 100 K = 1,NR
DO 100 I = 1,NL
IF (HED (I ,K,J) .LE. HEDSML) HED(I,K,J) = HEDSML
IF (HED(I,K,J) .GE. HEDLGE) HED(I,K,J) = HEDLGE
RINDEX = AL0G10(-1.*HED(I,K,J))*1.E01-1.1E01
INDEX = RINDEX
IF (HED( I ,K ,J) .LT. HEDT( INDEX+1)) INDEX = INDEX+1
IF (HED(I,K,J) .GT. HEDT(INDEX)) INDEX = INDEX-1
IF (INDEX .GE. 59) INDEX = 59
IF (INDEX .LE. 1) INDEX = 1
RM = RINDEX-INDEX
RATIO = AL0G10(HED( I, K,J)/HEDT( INDEX)) / AL 0G10 (HEDT( INDEX+1)/
+ HEDT(INDEX))
WCN(I,K,J) = WCNT( INDEX)+RATIO*(WCNT( INDEX+1)-WCNT( INDEX))
SWC(I ,K,J) = SWCT(INDEX)*1 ,E01**(RATIO*ALOG10(SWCT(INDEX+1)/
+ SWCT(INDEX)))
CON( I ,K ,J) = CONT( INDEX) *1 .E01**( RATI0*AL0G10(CONT( INDEX+1) /
+ CONT( INDEX)))
100 CONTINUE
RETURN
END
C
SUBROUTINE TBLYOL ( ILOOP )
C THIS SUBROUTINE USES AN INDEX ROUTING TECHNIQUE WITH
C TABULATED DATA TO FIND THE CORRESPONDING VALUES OF
C WATER CONTENT, SPECIFIC WATER CAPACITY AND HYDRAULIC
C CONDUCTIVITY FOR A GIVEN VALUE OF SOIL MATRIC POTENTIAL
161
C FOR YOLO LIGHT CLAY.
C - -
PARAMETER (IZ=50,IR=16,IB=IR*2+1,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
CHARACTER*64 MBFILE
REAL WCNT(120),CONT(120),SWCT(120),HEDT(120)
C0MM0N/PARA/WCN(IZ,IR,3),CON(IZ,IR,3) ,SWC(IZ,IR,3),HED(IZ,IR,3),
$ NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),ICONST,NR,DR(IR),R(IR),IL0W
COMMON/TABLES/WCNT,CONT,SWCT,HEDT,NDIM,MBFILE ,ISOIL
HEDSML = HEDT(49)
HEDLGE = HEDT(l)
DO 100 J = ILOOP ,3
DO 100 K = 1,NR
DO 100 I = 1,NL
IF (HED(I,K,J) .LE. HEDSML) HED(I ,K,J) = HEDSML
IF (HED(I ,K ,J) .GE. HEDLGE) HED(I,K,J) = HEDLGE
RINDEX = AL0G10(-1.*HED(I,K,J))*1.E01-1.1E01
INDEX = RINDEX
IF (HED(I,K,J) .LT. HEDT(INDEX+1)) INDEX = INDEX+1
IF (HED(I ,K,J) .GT. HEDT(INDEX)) INDEX = INDEX-1
IF (INDEX .GE. 49) INDEX = 49
IF (INDEX .LE. 1) INDEX = 1
RM = RINDEX-INDEX
RATIO = AL 0G10(HED(I,K,J)/HEDT(INDEX))/AL 0G10(HEDT(INDEX+1)/
+ HEDT(INDEX))
WCN( I ,K,J) = WCNT( INDEX) +RATI0* (WCNT( INDEX+1)-WCNT( INDEX))
SWC(I ,K,J) = SWCT(INDEX)*1,E01**(RATIO*ALOG10(SWCT(INDEX+1)/
+ SWCT(INDEX)))
CON(I,K,J) = C0NT( INDEX) *1 .E01**( RATI0*AL0G10(C0NT( INDEX+1) /
+ CONT(INDEX)))
100 CONTINUE
RETURN
END
C
SUBROUTINE LINEAR ( ILOOP )
PARAMETER (IZ=50,IR=16,IB=IR*2+1 ,MAX=IZ*IR)
IMPLICIT REAL (A-H.O-Z)
CHARACTER*64 MBFILE
REAL WCNT(120),CONT(120),SWCT(120),HEDT(120)
COMMON /TABLES/WCNT,CONT,SWCT,HEDT,NDIM,MBFILE,ISOIL
COMMON/PARA/WCN(IZ,IR,3),CON(IZ ,IR,3),SWC(IZ,IR,3),HED(IZ ,IR,3),
$ NL,DZ(IZ),DEPTH(IZ),QS(IR),S(IZ,IR),I CONST,NR,DR(IR),R(IR),ILOW
C
DO 100 K= ILOOP, 3
DO 100 J=1,NR
DO 100 I = 1,NL
IF(HED(I,J,K) .GT.O.O) HED(I,J,K) = 0.0
CALL INTPL( HED( I ,J ,K), CON( I ,J ,K) ,SWC(I ,J ,K) ,WCN(I ,J ,K))
100 CONTINUE
RETURN
END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE INTPL (W,X,Y,Z)
O C~) C“>
162
IMPLICIT REAL (A-H.O-Z)
THIS SUBROUTINE INTERPOLATES LINEARLY BETWEEN TWO SETS OF DATA
POINTS. DATA IS INPUT AS TABULAR FUNCTIONS
CHARACTER*64 MBFILE
REAL WTAB(120),XTAB(120),YTAB(120),ZTAB(120)
COMMON /TABLES/ZTAB,XTAB,YTAB,WTAB,NDIM,MBFILE,ISOIL
X= XTAB(1)
C Y=YTAB(1)
Z=ZTAB(1)
N=2
IF(ABS(W).LE.ABS(WTAB(1) )) GO TO 50
X=XTAB(NDIM)
C Y=YTAB(NDIM)
Z=ZTAB(NDIM)
N=NDIM
IF(ABS(W).GE.ABS(WTAB(NDIM) )) GO TO 50
C
DO 10 I=2,NDIM
N=I
IF(ABS(W).LE.ABS(WTAB(I))) GO TO 50
10 CONTINUE
50 X=XTAB(N-1)-ABS((W-WTAB(N-l))/(WTAB(N)-WTAB(N-1))*
+(XTAB(N)-XTAB(N-1)))
Y= YTAB(N— 1) -ABS((W-WTAB(N-l)) / (WTAB(N) -WTAB(N-1)) *
+(YTAB(N)-YTAB(N-l)))
Z=ZTAB(N-1)-ABS((W-WTAB(N-l))/(WTAB(N)-WTAB(N-l))*
+(ZTAB(N)-ZTAB(N-l)))
100 RETURN
END
REFERENCES
Amerman, C.R. 1969. Finite difference solutions of unsteady, two-
dimensional, partially saturated porous media flow. Volumes I and
II. Ph.D. Dissertation, Purdue University, W. Lafayette, Ind.
Amerman, C.R., and E.J. Monke. 1977. Soil water modeling II: On
sensitivity to finite difference grid spacing. Trans. Amer. Soc.
Agr. Eng. 20(3) :478-484,488.
Aribi, K., A.G. Smajstrla, F.S. Zazueta, and S.F. Shih. 1985.
Accuracies of neutron probe water content measurements for sandy
soils. Soil and Crop Sci. Soc. Fla. Proc. 44:44-49.
Armstrong, C.F., and T.V. Wilson. 1983. Computer model for moisture
distribution in stratified soils under a trickle source. Trans.
Amer. Soc. Agr. Eng. 26(6):1704-1709.
Bottcher, A.B., and L.W. Miller. 1982. Automatic tensiometer scanner
for rapid measurements. Trans. Amer. Soc. Agr. Eng. 15:1338-1342.
Brandt, A., E. Bresler, N. Diner, I. Ben-Asher, J. Heller, and D.
Goldberg. 1971. Infiltration from a trickle source: I.
mathematical models. Soil Sci. Soc. Amer. Proc. 35:675-682.
Bresler, E., J. Heller, N. Diner, I. Ben-Asher,A. Brandt, and D.
Goldberg. 1971. Experimental data and theoretical predictions.
Soil Sci. Soc. Amer. Proc. 35:683-689.
Clark, G.A. 1982. Simulation of the sensivity of water movement and
storage in sandy soils to irrigation application rate and depth.
Unpublished Master's Thesis. University of Florida, Gainesville,
FL.
Clark, G., and A.G. Smajstrla. 1983. Water distributions in soils as
influenced by irrigation depths and intensities. Soil and Crop Sci.
Soc. Fla. Proc. 42:157-165.
Denmead, O.T., and R.H. Shaw. 1962. Availability of soil water water
to plants as affected by soil moisture content and meteorological
conditions. Agron. J. 45:385-390.
Enfield, C.G., and C.V. Gillaspy. 1980. Pressure transducer for
remote data acquisition. Trans. Amer. Soc. Agr. Eng. 23(5): 1195-
1200.
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field
water use and crop yield. Halsted Press. New York.
163
164
Fitzsimmons, D.W., and N.C. Young. 1972. Tensiometer-pressure
transducer system for studying unsteady flow through soils.
Trans. Amer. Soc. Agr. Eng. 15:272-275.
Fletcher, J.E. 1939. A dielectric method of measuring soil moisture.
Soil Sci. Soc. Amer. Proc. 4:84-88.
Gardner, W.R. 1960. Dynamic aspects of water availability to plants.
Soil Sci. 89(2):63-73.
Gardner, W.R., and C.F. Ehling. 1962. Some observations on the
movement of water to plant roots. Agron. J. 54:453-456.
Gradmann, H. 1928. Untersuchungen uber die wasserverhaltnisse des
bodens ais grundlage des pflanzenwachstums. Jahrb. Wiss. Bot.,
69:1-100.
Hanks, R.J., and S.A. Bowers. 1962. Numerical solution of the moisture
flow equation for infiltration into layered soils. Soil Sci. Soc.
Amer. Proc. 26:530-534.
Harrison, D.S., A.G. Smajstrla, R.E. Choate, and G.W. Isaacs. 1983.
Irrigation in Florida agriculture in the '80s. Bulletin 196, Fla.
Coop. Ext. Serv., Univ. of Fla.
Haverkamp, R., M. Vauclin, J. Touma, P.J. Wierenga, and G. Vachaud.
1977. A compairison of numerical simulation models for one¬
dimensional infiltration. Soil Sci. Soc. Amer. Proc. 41:285-294.
Hiler, E.A., and S.I. Bhuiyan. 1971. Dynamic simulation of unsteady
flow of water in unsaturated soils and its application to
subirrigation system design. Technical Report TR-40. Texas Water
Resources Institute. College Station, TX.
Hi 11 el, D. 1980. Fundamentals of soil physics. Academic Press, New
York.
Hornberger, G.M., I. Remson, and A.A. Fungaroli. 1969. Numeric studies
of a composite soil moisture ground-water system. Water Resour.
Res. 5:797-802.
Long, F.L., and M.G. Hulk. 1980. An automated system for measuring
soil water potential gradients in a rhizotron soil profile. Soil
Sci. 129:305-310.
Marthaler, H.P., W. Vogelsanger, F. Richard, and P.J. Wierenga. 1983.
A pressure transducer for field tensiometers. Soil Sci Soc. Amer.
Proc. 47:624-627.
Molz, F.J., and J. Remson. 1970. Extraction term models of soil
moisture use of transpiring plants. Water Resour. Res. 6:1346-
1356.
Nimah, M.N., and R.J. Hanks. 1973a. Model for estimating soil water,
plant and atmospheric interrelations: I. description and
sensitivity. Soil Sci. Soc. Amer. Proc. 37:522-527.
165
Nimah, M.N., and R.J. Hanks. 1973b. Model for estimating soil water,
plant and atmospheric interrelations: II field test of model.
Soil Sci. Soc. Amer. Proc. 37:528-532.
Perrens, S.J., and K.K. Watson. 1977. Numerical analysis of two-
dimensional infiltration and redistribution. Water Resour. Res.
13(4):781-790.
Phene, C.J., 6.J. Hoffman, and S.L. Rawlins. 1971. Measuring soil
matric potential in situ by sensing heat dissipation within a
porous body: I. Theory and sensor construction. Soil Sci. Soc.
Amer. Proc. 35:27-33.
Phene, C.J., S.L. Rawlins, and G.J. Hoffman. 1971. Measuring soil
matric potential in situ by sensing heat dissipation within a
porous body: II. Experimental Results. Soil Sci. Soc. Amer.
Proc. 35:225-229.
Phillip, J.R. 1968. Steady infiltration from buried point sources and
spherical cavities. Water Resour. Res. 4:1039-1047.
Ritchie, J.T. 1973. Influence of soil water status and meteorological
conditions on evaporation from a corn canopy. Agron. J.
65:893-897.
Reddell, D.L., and D.K. Sunada. 1970. Numerical simulation of
dispersion in groundwater aquifers. Hydrology paper 41, Colorado
State University, Fort Collins, Colorado.
Rubin, J., and R. Steinhardt. 1963. Soil water relations during rain
infiltration: 1. Theory. Soil Sci. Soc. Amer. Proc. 27:246-251.
Rubin, J. 1967. Numerical method for analyzing hysteresis-affected,
post-infiltration redistribution of soil moisture. Soil Sci. Soc.
Amer. Proc. 31:13-20.
Rubin, J. 1968. Theoretical analysis of two-dimensional, transient
flow of water in unsaturated and partly unsaturated soils. Soil
Sci. Soc. Amer. Proc. 32:607-615.
Saxton, K.E., H.P. Johnson, and R.H. Shaw. 1974. Modeling
evapotranspiration and soil moisture. Trans. Amer. Soc. Agr. Eng.
17:673-677.
Slack, D.C., C.T. Haan, and L.G. Wells. 1977. Modeling soil water
movement into plant roots. Trans. Amer. Soc. Agr. Eng. 20:919-
927,933.
Smajstrla, A.G. 1973. Simulation of Miscible displacement in soils and
sensitivity to the dispersion coefficient. Master's Thesis. Texas
A & M University.
166
Smajstrla, A.G. 1982. Irrigation management for the conservation of
limited water resources. Office of Water Research and Technology
of the department of Interior Water Conservation Research Program
(WRC-80).
Smajstrla, A.G., 1985. A field lysimeter system for crop water use and
water stress studies in humid regions. Soil and Crop Sci. Soc. Fla.
Proc. 44:53-59.
Smajstrla, A.G., L.R. Parsons, F.S. Zazueta, G. Vellidis, and K. Aribi.
1986. Water use and growth of young citrus trees. ASAE Paper No.
86-2069, Amer. Soc. of Agr. Eng. St. Joseph, MI.
Taylor, G.S., and J.N. Luthin. 1969. Computer methods for transient
analysis of water-table aquifers. Water Resour. Res. 5:144-152.
Thomson, S.J., E.D. Threadgill, and J.R. Stansell. 1982.
Field test of a microprocessor based center pivot irrigation
controller. ASAE Paper No. 82-2535, Amer. Soc. of Agr. Eng.
St. Joseph, MI.
Thurnau, D.H. 1963. Algorithm 195, bandsolve. Communications of the
Association for Computing Machinery, 6(8):441.
Tollner, E.W., and F.J. Molz. 1983. Simulating plant water uptake in
moist, lighter textured soils. Trans. Amer. Soc. Agr. Eng.
26(1):87-91.
Van Bavel, C.H.M., G.B. Strik, and K.J. Burst. 1968. Hydraulic
properties of a clay loam soil and the field measurement of water
uptake by roots. I. Interpretation of water content and pressure
profiles. Soil Sci. Soc. Amer. Proc. 32:310-317.
Van den Honert, T.H. 1948. Water transport in plants as a catenary
process. Discuss. Faraday Soc. 3:146-153.
Wooding, R.A. 1968. Steady infiltration with a shallow circular pond.
Water Resour. Res. 4:1259-1273.
Zazueta, F.S., A.G. Smajstrla, and D.S. Harrison. 1984.
Microcomputer control of irrigation systems I: Hardware and
software considerations. Soil and Crop Sci. Soc. Fla. Proc.
43:123-129.
Zazueta, F.S., A.G. Smajstrla, and D.S. Harrison. 1985. A simple
numerical model for the prediction of soil-water movement from
trickle sources. Soil and Crop Sci. Soc. Fla. Proc. 44:72-76.
Zur, B., and J.W. Jones. 1981. A model for the water relations,
photosynthesis , and expansive growth of crops. Water Resour.
Res. 17:311-320.
BIOGRAPHICAL SKETCH
Kenneth Coy Stone was born June 25, 1959, in Tifton, Georgia. He
graduated from Tiftarea Academy, Chula, Georgia, in June 1977. He then
attended Abraham Baldwin Agricultural College in Tifton, Georgia, from
which he received an Associate in Science degree in mathematics in June
of 1979. In September of 1979 he began his undergraduate studies in
agricultural engineering at the University of Georgia. He graduated
from the University of Georgia with a Bachelor of Science in
Agricultural Engineering degree in June 1981. He then began his
graduate studies at the University of Georgia and received his Master
of Science (agricultural engineering) degree in March 1985. He enrolled
in the University of Florida in May 1983 to pursue the Ph.D. degree in
agricultural engineering..
167
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
ML i
D -
(l tIL
Allen G. Smajstrlaj
i Chairman
Associate Professor of
Agricultural Engineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
{il -
'ames W. Jones
fes sor of
Agricultural Engineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
Fadro R. Zaz
Assistant Pr
Agricultural
retai
Ifesápr of
ingineering
I certify that I have read this study and that in my opinion it
conforms to acceptable standards of scholarly presentation and is fully
adequate, in scope and quality, as a dissertation for the degree of
Doctor of Philosophy.
Assistant Professor of
Civil Engineering
This dissertation was submitted to the Graduate Faculty of the College
of Engineering and to the Graduate School and was accepted as partial
fulfillment of the requirements for the degree of Doctor of Philosophy.
May 1987
I'LsJLuJ' Q. ■
Dean, College of Engineering
Dean, Graduate School
UNIVERSITY OF FLORIDA
llllllllllllllllllllllllllllllllllllllllllllll
3 1262 08554 0135