MEASURED AND SIMULATED SOIL WATER
REDISTRIBUTION AND EXTRACTION PATTERNS OF
DRIPIRRIGATED TOMATOES ABOVE A SHALLOW WATER TABLE
By
GEORGE VELLIDIS
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
1989
Copyright 1989
by
George Vellidis
To my parents ..
who encouraged me to reach for the stars
ACKNOWLEDGEMENTS
I would like to express special appreciation and thanks to Dr. Alien G. Smajstrla for
his guidance, assistance, and encouragement throughout this research. His reliable support
and sound advice formed the basis of a rewarding and enjoyable experience at the University
of Florida. I would also like to thank the other members of my committee, Dr. D.Z. Haman,
Dr. F.S. Zazueta, Dr. G.A. Clark, and Dr. J.J. Delfino for their support.
Special thanks go to Dr. Ken Stone for his friendship and help and to Mr. Danny
Burch for his technical assistance and practical advice during the experimental phase of this
research.
Finally, I would like to thank my wife Tracey who helped me collect data, weigh
samples, irrigate tomatoes, repair shelters, and edit my writing. Her encouragement and her
beautiful smile made the problems I encountered much easier to overcome.
TABLE OF CONTENTS
tage
ACKNOWLEDGEMENTS .................. .................. .. iv
LIST OF TABLES ............. ........................ ............ vii
LIST OF FIGURES ............. .... ................... ............. ix
ABSTRACT ........... ....... ................... .... ............. xii
CHAPTERS
I INTRODUCTION ...................................... I
2 REVIEW OF LITERATURE .................. ........... 3
Soil W ater Potential ..................................... 3
Measurement of Soil Water .................................. 4
SoilMoisture Content Measurements .................. 4
Matric Potential Measurements ................. ....... 5
Automation of Matric Potential Measurements ........... 6
Soil Hydraulic Properties ................................ 8
SoilMoisture Retention Curve ................ ........ 8
Hydraulic Conductivity Curve ................. ....... 8
Water Flow in Unsaturated Soils ............................ 11
N umerical Solutions ................................ 14
Infiltration From a Line Source ...................... 15
Water Extraction by Plant Roots ........................... 17
Lysimetry ..... ............................... ........... 23
3 EXPERIMENTAL PROCEDURE ............................ 24
The Field Lysimeter System ............................... 24
Lysimeter Construction ............................ 24
The Data Acquisition System ........................ 31
Lysimeter System Operation ....................... 36
Soil Characteristics ........................................ 38
Bulk Density ...................................... 39
SoilMoisture Retention Curves ....................... 40
Saturated Hydraulic Conductivity ..................... 40
Hydraulic Conductivity Curve ......................... 41
4 EXPERIMENTAL RESULTS ................. ......... . 45
Plant G rowth ........................................... 45
M easurement Accuracy ........ ............................ 46
Water Budget Accuracy .............................. 46
Weekly, Monthly, and Seasonal Water Use Accuracy ........ 46
Crop Water Use ........................................ 47
Partitioning of Crop Water Use ....................... 47
Weekly, Monthly, and Seasonal Water Use ............... 49
Redistribution and Extraction Patterns ........................ 58
Evaporation from the Soil ................................. 66
5 NUMERICAL SOLUTION ................................ 69
OneDimensional Model Development ........................ 69
Boundary Conditions ............................... 73
Solution ........................................ 76
Tim e Step ....................................... 77
Updating of Soil Parameters .......................... 78
Mass Balance..................................... 78
OneDimensional Model Verification ........................ 79
TwoDimensional Model Development ....................... 84
Boundary Conditions ............................... 87
Solution ........................................ 91
Water Extraction by Plant Roots ....................... 92
TwoDimensional Model Verification ........................ 97
Simulation of the Lysimeter Study ........................... 102
Estimated Daily ET ................................ 105
The Root Distribution Term, RDT ..................... 106
Simulation Results ................................. 109
M odel Applications ...................................... 118
6 SUMMARY AND CONCLUSIONS .......................... 119
The Lysimeter System .................................... 119
Field Data ............................................. 120
Numerical Simulation .................................... 121
M odel Applications ...................................... 121
APPENDICES
A SOIL TEST REPORTS ................................... 122
B NUMERICAL MODEL PROGRAM LISTING ................. 134
C TABULATED SOIL PARAMETERS ........................ 181
REFERENCES.................................................... 192
BIOGRAPHICAL SKETCH .......................................... 197
vi
LIST OF TABLES
table pnag
3.1. Volumes of Irrigation Applications and Equivalent Depths Over the Lysimeter
Areas. ..................................................... 38
3.2. Characteristic Soil Properties. .................................... 39
3.3. Parameters for van Genuchten's Equation. ......................... 41
4.1. Canopy Dimensions (m) ..................... ........... ... 45
4.2. Monthly and Seasonal Water Use. ................................. 56
5.1. Measured and Modified Saturated Hydraulic Conductivity. .............. 110
A.I. Routine Test Report and Standard Fertilization Recommendations: Builder's
sand sample I ................................................ 123
A.2. Routine Test Report and Standard Fertilization Recommendations: Builder's
sand sample 2 ................................................ 124
A.3. Routine Test Report and Standard Fertilization Recommendations: Builder's
sand sample 3.................... ............................ 125
A.4. Routine Test Report and Standard Fertilization Recommendations: Astatula
f.s. sample I ........................ ....................... 126
A.5. Routine Test Report and Standard Fertilization Recommendations: Astatula
f.s. sample 2................................................. 127
A.5. Routine Test Report and Standard Fertilization Recommendations: Astatula
f.s. sample 3. .................................. .............. 128
A.7. Routine Test Report and Standard Fertilization Recommendations: Pomona
sand sample I.................. .............................. 129
A.8. Routine Test Report and Standard Fertilization Recommendations: Pomona
sand sample 2. ............................................... 130
A.9. Routine Test Report and Standard Fertilization Recommendations: Pomona
sand sample 3. .............................................. 131
A.10. Micronutrients, Organic Matter, Copper Toxicity Test, and Soluble Salts. .. . 133
C.I. Hydraulic Properties of Arredondo Fine Sand. ...................... 181
C.2. Hydraulic Properties of Rehevot Sand. ............................. 184
vii
C.3. Hydraulic Properties of Adelanto Loam. ............................ 186
C.4. Hydraulic Properties of Pachappa Loam. ............................ 188
C.5. Hydraulic Properties of Yolo Light Clay. ........................... 190
LIST OF FIGURES
figure page
3.1. Field layout of the lysimeter system. ............................... 25
3.2. Side view giving details of an individual lysimeter soil water status monitoring
system .................................................... 27
3.3. Float valve and sump used to control the water table level. .............. 29
3.4. End view of the float valve showing details of the mechanism used to mount
the valve onto a threaded stand. .................................. 29
3.5. Side view of a lysimeter showing the location of the 10 tensiometers. ....... 33
3.6. End view of a lysimeter showing the location of the 10 tensiometers. ........ 33
3.7. Calibration curves of output voltage versus pressure for transducer no. 4 .... 34
3.8. Fluctuations in atmospheric pressure for the week of June 1319. ......... 35
3.9. Arrangement of the three soils in the lysimeters. ...................... 36
3.10. Soilmoisture retention curves for the three soils used in this study. ....... 42
3.11. Fitted and experimental soilmoisture retention curves for the three soils used
in this study................................................. 43
3.12. Hydraulic conductivity curves for the three soils used in this study. ........ 44
4.1. Partitioning of crop water use in May. ............................ 48
4.2. Partitioning of crop water use in June. ............................ 49
4.3. Weekly crop water use in builder's sand, May, 1988. ................... 50
4.4. Weekly crop water use in builder's sand, June, 1988. ................... 51
4.5. Weekly crop water use in Astatula f.s., May, 1988. ................... 52
4.6. Weekly crop water use in Astatula f.s., June, 1988. ................... 53
4.7. Weekly crop water use in Pomona sand, May, 1988. .................... 54
4.8. Weekly crop water use in Pomona sand, June, 1988. ................... 55
4.9. Cumulative water use for the plants growing in builder's sand ............. 56
4.10. Cumulative water use for the plants growing in Astatula f.s. ............ .57
4.11. Cumulative water use for the plants growing in Pomona sand. ............ 57
4.12. Tensiometer readings versus time on a daily scale. ................... .. 59
4.13. Tensiometer readings versus time on an hourly scale. ................... 59
4.14. Tensiometer readings versus time on a daily scale for the week of June 6. ... 61
4.15. Side view of a lysimeter showing the bed, the water table, and the tensiometer
locations.................................................... 61
4.16. Equipotential lines in a Pomona sand profile at irrigation (kPa). ........... 63
4.17. Equipotential lines in a Pomona sand profile during redistribution (kPa) .... 63
4.18. Equipotential lines in a Pomona sand profile during redistribution (kPa) .... 64
4.19. Equipotential lines in a Pomona sand profile following redistribution (kPa). 64
4.20. Equipotential lines in a Pomona sand profile during extraction (kPa). ........ 65
4.21. Equipotential lines in a Pomona sand profile during extraction (kPa). ........ 65
4.22. Comparison of Penman ET and observed bare soil evaporation from the
lysimeters ................................................... 67
4.23. Comparison of Penman ET and estimated daily bare soil evaporation from the
lysimeters during the final two months of the growing season. ............ 68
5.1. The finite difference grid system used for the onedimensional model.
Adapted from Stone (1987), p.38. ................................. 70
5.2. Simulated, analytical, and experimental results of steadystate evaporation at
a rate of 0.033 mm/hr from a water table in Yolo light clay. ............. 80
5.3. Simulated results of soil water content profiles for infiltration into Rehevot
sand under constant rain intensity of 12.7 mm/hr. ...................... 81
5.4. Simulated results of soil water content profiles for infiltration into Rehevot
sand under constant rain intensity of 47 mm/hr. ..................... 81
5.5. Simulated and analytical results of steadystate infiltration above a shallow
water table at rates of 0.0 and 1.5 mm/hr. ............................ 83
5.6. Simulated and analytical results of steadystate infiltration above a shallow
water table with a constant extraction rate of 0.42 mm/hr. ................ 83
5.7. Schematic diagram of the finitedifference grid system for the two
dimensional model of water movement and extraction. ................. 85
5.8. Effect of the relative available soil water and four ET rates on the predicted
soil water extraction rate ........................................ 95
5.9. Simulated distribution of the daily ET over a 24hr period. .............. 96
5.10. Twodimensional model simulated results of soil water content profiles for
infiltration into Rehevot sand under constant rain intensity of 12.7 mm/hr. .. 98
5.11. Twodimensional model simulated results of soil water content profiles for
infiltration into Rehevot sand under constant rain intensity of 47 mm/hr. ... 98
5.12. Lines of constant pressure head, 0 (mm), for parallel buried line sources
1220 mm apart ............................................... 100
5.13. Lines of constant pressure head, b (mm), for a buried line source with a water
table 300 mm below the soil surface................................ 100
5.14. The vertical water content distribution at (0,z) for a line source discharge of
11442 mm2/hr after 0.40 hrs. .................................... 103
5.15. The horizontal water content distribution at (x,210) for a line source discharge
of 11442 mm'/hr after 0.40 hrs................................... 103
5.16. Finite difference mesh for the lysimeter problem. ................... .. 104
5.17. Comparison of Penman ET and estimated daily ET from the lysimeters during
the final two months of the growing season. ........................ 105
5.18. The number of cells contained by each soil compartment. ............... 107
5.19. Distribution of extraction in Astatula f.s. for the week of May 2 ........... 108
5.20. Distribution of extraction in Pomona sand for the week of June 6. ......... 108
5.21. Comparison of simulated and observed pressure heads at a depth of 100 mm
in Astatula f.s. for the week of May 2. ............................. 111
5.22. Comparison of simulated and observed pressure heads at a depth of 200 mm
in Astatula f.s. for the week of May 2 ............................. 112
5.23. Comparison of simulated and observed pressure heads at a depth of 300 mm
in Astatula f.s. for the week of May 2. ............................ 113
5.24. Comparison of simulated and observed pressure heads at a depth of 100 mm
in Pomona sand for the week of June 6. ........................ 114
5.25. Comparison of simulated and observed pressure heads at a depth of 200 mm
in Pomona sand for the week of June 6. ............................ 115
5.26. Comparison of simulated and observed pressure heads at a depth of 300 mm
in Pomona sand for the week of June 6. .......................... 116
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
MEASURED AND SIMULATED SOIL WATER
REDISTRIBUTION AND EXTRACTION PATTERNS OF
DRIPIRRIGATED TOMATOES ABOVE A SHALLOW WATER TABLE
By
GEORGE VELLIDIS
May, 1989
Chairman: Allen G. Smajstrla
Major Department: Agricultural Engineering
In Florida, the greatest use of fresh water resources is for supplemental irrigation of
agricultural crops. Because of increasing competition for the consumptive use of this water,
there has been a rapid increase in the use of micro irrigation, and in particular, line source
drip irrigation systems which are increasingly being used to irrigate high cash value crops
such as fruits, vegetables, and ornamentals.
Throughout Florida, vegetables are often grown on mulched beds in soils that
commonly require irrigation even if they are subject to naturally occurring high water tables.
The contribution of high water tables to soil water extracted by agricultural crops had not
been accurately determined, especially in combination with line source drip irrigation
systems. This soilplantwater relationship problem lent itself well to a study conducted in
lysimeters.
A set of six rainsheltered field lysimeters was developed to simulate these conditions.
Each drainage lysimeter contained a water table regulation device and was irrigated with a
line source drip irrigation system. Soil pressure potential data were collected continuously
using a microcomputerbased tensiometerpressure transducer data acquisition system.
Experiments were conducted in the lysimeters to measure crop water use and
redistribution and extraction patterns of a tomato crop. The design of the lysimeters
permitted crop water use to be partitioned into irrigation and water table recharge
components. Results indicated that a significant portion of daily plant water use was
provided by the water table. Redistribution and extraction patterns calculated from the soil
water potential measurements showed that very little change in soil water status took place
at horizontal distances of more than 200 mm from the irrigation lateral and at depths greater
than 300 mm.
A twodimensional finite difference model was developed to simulate water
infiltration, redistribution, and extraction under the above conditions. The simulated results
were in excellent agreement with the field observations. The model can be used to
investigate many different irrigation strategies without the cost normally associated with field
experiments.
CHAPTER 1
INTRODUCTION
The greatest use of Florida's fresh water resources is for supplemental irrigation of
agricultural crops. Harrison et al. (1983) reported that 41% of the state's total fresh water use
was for irrigation. In 1987, 728,000 ha were irrigated in Florida by all types of irrigation
systems (Goldstein, 1988).
There are three main methods used for supplemental irrigation in Florida: sprinkler,
seepage, and micro irrigation. In the past 10 years there has been a rapid increase in the use
of micro irrigation in Florida. This is due to energy and water savings offered by these
systems in a time of increasing fuel costs, the ability to apply nutrients and chemicals directly
to the crop, and increasing competition for the consumptive use of available fresh water. Of
special interest are line source drip irrigation systems which are increasingly being used to
irrigate high cash value crops such as fruits, vegetables, and ornamentals.
Line source systems consist of irrigation laterals with closelyspaced discharge points
or emitters. These differ from other micro irrigation systems in that they generate a wetted
band of soil along the lateral, as opposed to a series of nonoverlapping wetted bulbs. When
well designed and managed, and in absence of rainfall, they provide adequate water for crop
growth and minimize losses from percolation below the root zone and evaporation from the
soil surface.
Florida soils commonly require irrigation even if they are subject to naturally
occurring high water tables (0.51.5 m below the bed surface). The contribution of high
water tables to soil water extracted by agricultural crops has not been accurately determined,
especially in combination with line source drip irrigation systems. This soilplantwater
relationship problem lends itself well to studies conducted in lysimeters if the lysimeter
system is designed to adequately approximate the physical system (Chow, 1964).
A sound knowledge of the dynamics of water movement from line sources into the
soil, and of the extraction patterns observed under agricultural crops is needed to minimize
irrigation input while optimizing production returns. The dynamics of interactions must be
studied mathematically using numerical methods becauseanalytical solutions are applicable
to only a few specific cases due to the nonlinear nature of the governing equations. Many
researchers have developed numerical models to use in the study of soil water movement
(Smajstrla, 1982; and Feddes et al., 1978). These models allow many different irrigation
strategies to be investigated without the cost normally associated with field experiments
(Stone, 1987). These models, however, do not address the unique conditions often
encountered in Florida.
The purpose of this research was to obtain a clear understanding of the dynamic soil
plantwater relationships that occur when agricultural crops are grown on sandy soils,
irrigated by line source drip systems, and subjected to naturally occurring high water tables.
To achieve this goal, three objectives were established:
1. To develop and test a set of field lysimeters in which irrigation studies could
be conducted under high water table conditions.
2. To use realtime monitoring of soil water potentials to record the dynamics
of soil water redistribution and extraction under line source micro irrigation
systems and high water table conditions.
3. To develop a mathematical model that describes the redistribution and
extraction of soil water governed by the above conditions and formulate a
numerical procedure for evaluating the model.
CHAPTER 2
REVIEW OF LITERATURE
Soil Water Potential
In the soil, water moves from points where it has a high energy status to points where
it has a lower one. It is constantly in pursuit of that elusive state known as equilibrium
(Hillel, 1980). Therefore, soil water moves constantly in the direction of decreasing potential
energy. The energy status of water in soil is called the water potential T and is composed
of several components
T= + Vg + ( osm+ gas) (2.1)
where
= matric potential, resulting from capillary and adsorptive forces due to the soil
matrix,
og = gravitational potential, resulting from the gravitational force,
osm = osmotic potential, resulting from osmotic forces, and
gas = pneumatic potential, resulting from changes in ambient gas pressure.
In studies of soil moisture flow, the osmotic and pneumatic potentials can usually be
neglected. The influence of 0oa is low because, in general, solute concentrations in the soil
solution will not significantly affect liquid flow. Pneumatic pressure can be disregarded
because atmospheric pressure remains nearly constant and gas pressures in natural soil
generally do not differ from the atmospheric pressure. Equation (2.1) then reduces to
= + og (2.2)
Soil water potential can be expressed in terms of mass, volume, or weight, and it is
most often expressed as energy per unit weight of soil water. Soil water potential then has
units of length, i.e. mm, which can easily be translated to hydrostatic pressure. For example,
100 mm of water is equivalent to about 1 kPa or 0.01 bar.
The matric potential (f) in unsaturated soil is negative, because work is needed to
withdraw water against the soil matric forces. At the phreatic surface, V=0 mm, while it is
positive below the water table.
The gravitational potential (0 ) at each point is determined by the height of the point
relative to some arbitrary reference level, z. If the origin of z is taken at the soil surface
with z positive in the downward direction, 0==z mm anywhere below the soil surface.
When dealing only with the sum of matric and gravitational potential expressed in
units of length, the sum is referred to as the hydraulic head, H. Equation (2.2) becomes
H z (mm) (2.3)
with 0 now called the soil moisture pressure head and z the gravitational head. It is
convenient to refer to negative pressure head (0) as a positive suction or tension (h). Thus
h = 0 (2.4)
and the value of h ranges from 0 to 108 mm (Feddes et al., 1978).
Measurement of Soil Water
There are two general reasons for measuring soil water. One is to determine the
moisture content of the soil (0), and the other to determine the magnitude of the matric
potential, or pressure head (t), in the soil.
SoilMoisture Content Measurements
There are direct or indirect methods to measure 0. The gravimetric method involves
removing a sample by augering into the soil and then determining the moist and dry weights
(Hillel, 1980).
Neutron scattering allows the nondestructive measurement of soilmoisture content.
Although the method is efficient and reliable, it has low spatial resolution. Gamma ray
absorption, on the other hand, allows measurement in very thin layers of soil but is a rather
difficult and erratic procedure.
Another approach involves the use of the electrical resistance of a soil volume to
measure its moisture content. Colman and Hendrix (1949), among others, used electrical
resistance blocks placed in the soil and left to equilibrate with soil moisture to produce a
calibration curve of block resistance versus soil wetness. They found, however, that porous
blocks placed in the soil equilibrate with matric potential rather than soil moisture.
Except for the gravimetric method, the techniques discussed produce indirect
measurements of soilmoisture content and do not lend themselves to automation. Overall,
the measurement of soil wetness does not provide a sufficient description of the state of soil
water (Hillel, 1980). To obtain such a description, matric potential must be measured.
Matric Potential Measurements
To measure matric potential in the field, an instrument known as a tensiometer is
used. The tensiometer consists of a closed tube with a porous ceramic cup on the end which
is inserted into the soil. The cup/tube assembly is filled with deaerated water, sealed, and
allowed to equilibrate with the soil water. Initially, moist soil at atmospheric pressure is
placed in contact with pure water across the ceramic. Solutes then diffuse from the soil
across the ceramic to the tensiometer solution until solute equilibrium is established. At that
time, the solution in the instrument is called the equilibrium dialyzate (Taylor and Ashcroft,
1972). This effectively eliminates any osmotic potential component from consideration. Soil
water, being generally at subatmospheric pressure, exercises a suction which draws out a
certain amount of dialyzate from the rigid and airtight tensiometer, thus causing a drop in
its hydrostatic pressure (Hillel, 1980). A pressure sensing device, usually a vacuum gage, a
mercury manometer, or an electrical transducer, is used to measure the pressure in the
tensiometer tube.
In a tensiometer, the matric potential is not measured directly. Instead, the negative
pressure created within the tensiometer by the matric suction is measured. When the pressure
head is divided by the density of liquid water, the pressure potential of the dialyzate in terms
of energy per mass is obtained. Although the matric potential of the soil is not measured
directly, it is assumed to be equal to the pressure potential of the dialyzate, which is directly
determined (Taylor and Ashcroft, 1972). Tensiometers left in the soil for a long period of
time tend to follow the changes in the matric potential. This potential is also the hydraulic
potential that a plant must exert to extract water from the soil.
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 changes in the soil
matric suction. Lags are also caused by the volume of water needed to be moved through the
cup to register on the measuring device. This lag time can be minimized by the use of a
nulltype device or of a transducertype pressure sensor with rigid tubing, so that practically
no flow of water is required as the tensiometer adjusts (Hillel, 1980).
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 crops on sandy soils
because 75% or more of the amount of soil water used by plants is depleted within the
tensiometer range.
Automation of Matric Potential Measurements
Early methods of measuring tensiometer pressure heads used mercury manometers.
Later, mechanical vacuum gauges were used. Both methods functioned well for manual
applications but neither was readily automated (Stone, 1987). Recent interest in better
understanding the dynamics of soil water movement has spurred the development of
numerical models that simulate water flow through the soil. The large amounts of data
required to validate and calibrate these models have resulted in the need for an automated
system of recording tensiometer readings on a continuous realtime basis.
Most recent efforts have focused on using a transducertype pressure sensor to read
the pressure head. Marthaler et al. (1983) used a pressure transducer to read individual
tensiometers. The upper end of the tensiometer was closed off with septum stoppers to
provide an airtight seal. A needle connected to a pressure transducer was manually inserted
through the septum into a pocket of entrapped air in the tensiometer and the pressure
transducer output was obtained immediately. However, since the transducer had to be moved
manually from tensiometer to tensiometer, this was not an automated system.
Thomson et al. (1982) used individual pressure transducers on tensiometers with all
air purged from the system to monitor matric potentials. Thus, they were able to avoid lag
times associated with air pockets in the tensiometers. They read the pressure transducers by
electronically switching from one to another.
Stone et al. (1985) developed a lowcost microcomputerbased data acquisition system
for continuous pressure head measurements. Their system consisted of tensiometermounted
pressure transducers, a multiplexing system, an analogtodigital converter, and a portable
microcomputer. The system components were calibrated independently and the assembled
data acquisition system was evaluated under laboratory and field conditions. They achieved
excellent agreement between pressure heads read with the data acquisition system and those
read manually using mercury manometers. Agreement was generally within 0.47 kPa. The
system was tested under various field conditions and proved to be reliable and accurate.
Although the system developed by Stone et al. (1985) was inexpensive and reliable,
it involved developing most of the system components. In the time since then, numerous data
acquisition and control devices have become commercially available. These devices receive
an analog signal from the instruments, convert the signal to digital form, and pass it on to
a microcomputer. In addition, most of the devices have up to 32 channels available for data
input. This enhanced ability allows a powerful microcomputer based data acquisition system
for continuous soil water pressure measurement to be assembled from commercially available
instruments.
Soil Hydraulic Properties
Solution of unsaturated flow problems requires the predetermination of the soil
hydraulic properties, namely, the relationship between the pressure head, i, and the soil
moisture content, 0, and the dependence of the hydraulic conductivity, K, on 0 or 0. The
most important single factor limiting the successful application of unsaturated flow theory
to actual field problems is the lack of information on the parameters entering the governing
flow equations (van Genuchten, 1980).
SoilMoisture Retention Curve
The amount of water retained in the unsaturated zone at equilibrium is a function of
the sizes and volumes of the waterfilled pores. It is also, therefore, a function of the matric
forces. This function can be determined experimentally, and is represented graphically by
a curve known as the soilmoisture retention curve (Childs, 1940).
The soilmoisture retention curve is a plot of pressure head (V) versus soilmoisture
content (6) and is obtained by desorption or sorption. The equilibrium moisture content at
a given 0 is greater in desorption (drying) than in sorption (wetting) (Hillel, 1980). This
phenomenon is referred to as hysteresis. Poulovassilis (1962) presented a rigorous
mathematical description of hysteresis. In many studies hysteresis is not considered since a
soil's spatial variability is unknown and often exceeds the influence of hysteresis (Feddes et
al., 1978). Desorption curves are also much easier to develop and are therefore typically used
to describe the soilmoisture retention curve.
The desorption curve is usually determined by using a gas pressure extractor which
provides a reliable means of removing soil moisture from soil samples, under controlled
conditions, throughout the whole plant growth range.
Hydraulic Conductivity Curve
Hydraulic conductivity is defined as the ratio of the flux, q, to the hydraulic gradient,
or the slope of the flux versus gradient curve (Hillel, 1980). For saturated flow, the total soil
9
pore space is available for water flow. With unsaturated flow, however, part of the pores are
filled with air. Therefore, the unsaturated hydraulic conductivity, K, is smaller than the
saturated hydraulic conductivity, K,. So for unsaturated soils, K is not constant but depends
on the soil moisture content 0 or the pressure head 0 (Feddes et al., 1978).
K = f(0) or K = f(V) (2.5)
Direct measurement of unsaturated hydraulic conductivity is desirable but difficult
to obtain because experimental determination is timeconsuming and costly, relationships
between K, 0, and 6 are hysteretic and dependent upon previous wetting and drying histories,
and soil variability is such that an enormous amount of data would be required to represent
the hydraulic conductivity accurately (Mualem, 1986). To replace the missing information,
several models have been developed to compute hydraulic conductivity. Mualem (1986)
presented a review of several models.
One approach that has been successful involves calculating the unsaturated
conductivity from the more easily measured soilmoisture retention curve. The Millington
Quirk method (Millington and Quirk, 1971) has been applied with good results by several
researchers including Green and Corey (1971). The method's disadvantage is that it produces
output in tabular form which is difficult to use when dealing with nonhomogeneous, multi
dimensional, unsaturated flow problems (van Genuchten, 1980).
Brooks and Corey (1966) developed a closedform analytical expression for predicting
unsaturated hydraulic conductivity whose predictions compared favorably with experimental
results. However, a discontinuity which prevented rapid convergence in numerical
simulations of unsaturated flow problems was encountered in their soilmoisture retention
and hydraulic conductivity curves.
Mualem (1976) derived a simple integral formula from which closed form analytical
solutions for predicting the hydraulic conductivity were derived, provided suitable equations
for the soilmoisture retention curves were known. Van Genuchten (1980) derived such
closedform analytical solutions. The conductivity models he developed contained three
independent parameters which were obtained by matching a proposed soilmoisture retention
function to experimental data. Agreement was excellent between observed and predicted
conductivity curves for lighter soils such as sands and loams but higher conductivity values
were seriously underpredicted for a clay soil.
The soilmoisture content as a function of the pressure head is given by
van Genuchten (1980) as
(9, 0)
S= ,O + (2.6)
I1 + (o lI )"nJm
where
0, = saturated soilmoisture content,
Or = residual soilmoisture content,
10 = absolute value of the pressure head,
a, n = parameters to be determined, and
m = 1 I/n. (2.7)
Equation (2.6) contains four independent variables (0,, 0r, a, and n), which must be
obtained from observed soilmoisture retention data. The most available of the four is 0,,
which is easily determined experimentally. The residual water content, 0r, may be measured
experimentally by determining the soilmoisture content on very dry soil. It is defined as
the low end soilmoisture content for which the gradient (dB/dt) becomes zero, excluding
the region near 0,, which also has a zero gradient (van Genuchten, 1980). The remaining
parameters, a and n, are determined by differentiating equation (2.6) to obtain
dO am(0, 0,) 1 (2.8
 \ (2.8)
d ]m +(a J L +(a )
where the right hand side is expressed in terms of 0. Further substitutions solve for a, n, and
m. A step by step procedure was provided by van Genuchten (1980).
II
The hydraulic conductivity was expressed in terms of pressure head by the following
expression
( (alt I I)"' + (al I )~ 2
K(t) = K, (2.9)
[I +(all )"]m/2
where all variables are as defined for equation (2.6).
Water Flow in Unsaturated Soils
The general equations governing unsaturated flow in porous media are Darcy's law
and the continuity of flow equation. Richards (1931) linked the two and presented a
combined flow equation which he described as the general partial differential equation of soil
water transfer.
Darcy's law for the flow of water in unsaturated porous media (Darcy, 1856),
including soils, can be expressed as
q = K(b)VH (2.10)
where
q = density flux or volume of water passing through a unit crosssectional area per
unit time,
K = unsaturated hydraulic conductivity of the soil as a function of 0,
b = soil moisture pressure head,
V = differential operator, and
H = hydraulic head.
The continuity of flow equation according to Freeze and Cherry (1979) can be
expressed as
80
Vq S (2.11)
at
where
0 = soilmoisture content,
t =time, and
S = a sink or source term.
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.
Substitution of equation (2.10) into the continuity equation yields the combined flow
equation or the Richards equation which can be expressed as
80
 V[ K(O) VH] S() (2.12)
Bt
Full expansion of equation (2.12) in two dimensions yields the flow equation in terms
of hydraulic head
= K(4) + K() I S(4) (2.13)
at ax L ax 8z 8z
where x is the horizontal dimension and z is the vertical dimension with z positive in the
downward direction.
Substitution of equation (2.3) into equation (2.13) yields the pressure head form of
the flow equation
=  K() + K() S() (2.14)
at ax ax az zL
Differentiating the internal terms of equation (2.14) with respect to x and z
accordingly yields
86 8 a 8( ) 8 [ f8(4) 11
= K() + K() S(4) (2.15)
at ax L x J a I az
Childs (1957) used a different approach to obtain the same general flow equation.
Equation (2.15) is complicated by the appearance of two dependent variables 0 and
V. Rubin (1966) discussed three ways of reducing the equation to a single dependent variable,
among which was the Obased approach which was selected for this study. This technique
involves modifying equation (2.15) by introducing the specific water capacity, C, defined by
dO
C = (2.16)
do
which was discussed in detail by Klute (1965). Because hysteresis is not included, is
considered a singlevalue function of 0 and a normal instead of a partial derivative notation
is used in equation (2.16) (Feddes et al., 1978).
Using the chain rule of calculus
80 dO a8
(2.17)
at d 8at
and substituting equation (2.16) into (2.17) yields
8a 84
= C (2.18)
8t at
which in turn is substituted into equation (2.15) to produce the final form of the general flow
equation, commonly known as the two dimensional form of the Richards equation.
84 a [ 8(a) 8(0) )I
C() K() + K(') 1 S) (2.19)
a1 8x L x zL tz
Written in this form, the flow equation is easier to solve and provides the basis for
predicting soil water movement in saturated and unsaturated soils. The equation is a
nonlinear, parabolic, second order, partial differential equation (PDE). It is nonlinear because
the parameters K(4), C(Q), and S(O) depend on the solution of O(x,z,t). Nonlinearity causes
problems in the equation's solution and analytical solutions are complex and known for special
cases only (e.g., Gilley and Allred, 1974; Warrick et al., 1979). The nonlinearity of equation
(2.19) and the typical variable boundary conditions associated with it have dictated the use
of numerical methods to solve practical field problems of soilplantwater relationships
(Stone, 1987).
Numerical Solutions
The finite difference method is the numerical technique most commonly used to solve
the Richards equation. Finite differences are expedient because solutions are obtained by
working directly with the governing differential equations, and they apply to specific flow
geometries with specific mesh spacings. Therefore, finite differences lend themselves well
to the solution of soilplantwater relationship problems in regularly shaped lysimeters.
Finite difference formulations may be solved explicitly or implicitly. Although the
implicit approach is more complicated, it is preferable because it offers better stability and
convergence. It also permits relatively larger time steps, thus allowing longer simulation
periods.
The Richards equation has been successfully solved for onedimensional flow by
many researchers. Haverkamp et al. (1977) reviewed six numerical models of one
dimensional infiltration solved using both the water content based equation and the pressure
head based equation. Each model employed different discretization techniques to
approximate the nonlinear infiltration equation. They found that implicit models which
solved the pressure head based infiltration equation had the widest range of applicability for
predicting water movement in saturated or unsaturated soil.
Feddes et al. (1978) developed a onedimensional implicit model that simulated
infiltration and redistribution under high water table conditions. Plant growth and root
extraction were also simulated. Field lysimeters were used to collect data for model
calibration and special attention was paid to the influence of the groundwater table on the
amount of water available for transpiration of red cabbage on sticky clay.
15
The most popular implicit procedure used in solving the twodimensional form of the
Richards equation has been the alternating direction implicit (ADI) method. Rubin (1968)
developed a twodimensional numerical model of transient water flow in unsaturated and
partly unsaturated soils which utilized an ADI scheme. He studied horizontal infiltration and
ditch drainage with the numerical model.
Hornberger et al. (1969) developed a model of transient water movement in a
composite soil moisture ground water system. The model was used to study twodimensional
flow in response to a falling water table and used the GaussSidel iterative method to solve
the flow equations.
Amerman (1976) developed a computational scheme for solving the steadystate, two
dimensional form of the Richards equation using the successive overrelaxation (SOR) finite
difference method. The model was developed using subsectioning which allowed application
of the model, without internal modification, to a wide range of geometric shapes, hydraulic
boundary conditions, and soil distributions.
Perrens and Watson (1977) studied a twodimensional infiltrationredistribution
sequence in which the surface flux was spatially nonuniform. They found that an ADI
method described the flow process satisfactorily. Their program incorporated an
interpolationtype hysteresis model and produced twodimensional profile distributions of
water content and pressure head.
Brandt et al. (1971) considered a plane flow and a cylindrical flow model in a study
of twodimensional transient infiltration from a trickle source. The water flow equation was
solved numerically by an approach that combined the noniterative ADI procedure with
Newton's iterative method. The numerical results compared favorably with data from the
literature and indicated that the method was reliable and could be used with confidence.
Infiltration From a Line Source
A special twodimensional water flow problem is infiltration from a line source.
Researchers have presented several analytical and numerical solutions to this problem. Raats
(1970) was one of the earliest investigators to study steady infiltration from an array of
equally spaced line sources or furrows at the surface of a semiinfinite soil profile. He used
the steadystate assumption and the Kirchhoff transformation to obtain a straightforward
analytical solution. He produced explicit expressions for the stream function, the flux, the
matric flux potential, the pressure head, and the total head.
In a similar steadystate study, Zachmann and Thomas (1973) developed equations
which can be used in the design of subsurface irrigation systems that lie in a horizontal plane
and are parallel and equally spaced. The analytical solution described the flow from the line
sources in the presence of uniform infiltration or evaporation at the soil surface.
Warrick and Lomen (1977) developed a twodimensional analytical solution of the
linearized flow equation for a buried line source with an underlying constant potential
surface. Their geometry was appropriate for describing flow from a trickle or subsurface
irrigation system of porous pipe, or emitters, spaced closely in a line above a shallow water
table or in a lysimeter. They examined only the steadystate case and attained linearization
by assuming that the unsaturated hydraulic conductivity was exponentially related to the
pressure head. Warrick et al. (1979) further enhanced the model by adding a onedimensional
extraction term. Other analytical solutions of infiltration from a line source were provided
by Gilley and Allred (1974), Lomen and Warrick (1974), and Thomas et al. (1974).
Fewer studies have provided numerical solutions to the line source problem. Oron
(1981) linearized the steadystate nonlinear equation and then solved it numerically. His
model included a continuous sink function which simulated water uptake by roots. Oron's
results satisfactorily simulated the flow characteristics in the soil and were validated by
experimental data.
Ragab et al. (1984) used a simulation model based on the matric flux potential concept
to describe twodimensional infiltration into sand from a line source. The model showed
considerable stability with varying grid sizes and flow rates. Their simulated wetting patterns
indicated close agreement with published experimental data, especially under higher flow
rates.
An extensive literature review has shown that many researchers have investigated
infiltration from line sources. Most however, limited their studies to the solution of the
linearized, steadystate form of the flow equation. Those that solved for transient, two
dimensional infiltration from line sources did not include an extraction term to simulate
water uptake by plants. Furthermore, the line source problem has not been investigated
under high water table conditions. A general model simulating transient infiltration and
extraction under line sources is lacking.
Water Extraction by Plant Roots
Hamon (1966) calculated that 70% of the 750 mm average annual precipitation of the
contiguous United States is lost to evaporation and transpiration. Penman (1970) estimated
that 2000 Mg of water must be extracted from the soil in order to produce 20 Mg of fresh
weight of a crop. These observations indicate that the soilplantatmosphere pathway is a
major component of the subsurface hydrologic system. Any serious attempt to model the
pathway requires an understanding of the dynamics of soil moisture depletion by plants.
Most models of soil extraction by plant roots use the Ohm's law analogy which was
originally introduced by Gradmann (1928) and enhanced by van den Honert (1948). The
analogy suggests that the driving force for water uptake is the effective potential gradient
between water in the plant and water in the soil. The resistance term is a summation of the
hydraulic resistances of the soil as well as the root cortex and xylem along the flow path
(Belmans et al., 1979). In simulating water extraction by plant roots, two methods have been
used: the microscopic and the macroscopic approach.
A microscopic study examines the flow process in the vicinity of a single root. The
root may be considered as a cylindrical sink of infinite length, uniform radius, and uniform
water absorbing properties. The soil moisture flow equation is written in cylindrical
coordinates and solved with appropriate boundary conditions at the root surface and at some
distance rm, from the root (Molz and Remson, 1970). Moisture flux, moisture content, or
a combination of both, may be specified as a boundary condition at the root surface. Gardner
(1960; 1964) developed one of the earliest detailed quantitative models of this type. Other
microscopic studies were conducted by Molz et al. (1968), Newman (1969), Whisler et al.
(1970), and Steinhardt et al. (1982).
As discussed by Molz and Remson (1970), microscopic models are impractical because
flow to each individual rootlet of a complete root system must be considered. Furthermore,
the entire root system is not uniformly permeable to water and its geometry is time
dependent. As a result, most extraction functions have been developed using a macroscopic
approach (Molz, 1981).
In the macroscopic approach, the flow to individual roots is disregarded and the entire
root system is assumed to extract moisture from each differential volume, or cell, of the root
zone at some rate S. At a given cell, the quantity S is interpreted as an extraction rate that
may be a function of space, time, moisture content, or a combination of these variables
(Slack, 1975). The total extraction rate is the summation of the water flow rates from each
soil cell. In the soil moisture flow equations, the extraction rate is represented as a sink or
a negative source term, S, and is simply added to the continuity equation (see equation 2.11).
Boundary conditions are specific at boundaries of the composite soilplant system such as the
soil surface and the water table, and they do not involve the actual roots (Molz and Remson,
1970).
Gardner and Ehlig (1962) used the Ohm's law analogy to present a macroscopic
extraction function of the form
S = pn (2.20)
Ra01 + Rroot.
where
S = the extraction function or the volume of water taken up by the roots per unit
bulk volume of soil in unit time,
fsili ,plant = water potentials in the soil and plant respectively, and
Roil, Ro,,, = the corresponding hydraulic resistances.
They conducted greenhouse experiments and found that the impedance to water movement
in the soil is greater than the impedance to water movement into plant roots. Their
experimental results were consistent with the uptake equation.
Molz and Remson (1970) developed a onedimensional model which used a
macroscopic extraction term to describe moisture removal from soil by the roots of
transpiring plants. The extraction rate was given by
R(z)D(0)
S(z,O) = T (2.21)
JR(z)D(0)dz
where
T = transpiration rate per unit area,
R(z) = effective root density,
D(O) = diffusivity, and
v = root depth.
A numerical procedure based on the DouglasJones predictorcorrector method was used to
solve the model. Solutions compared with experimental data indicated that the extraction
term was computationally and physically feasible.
Belmans et al. (1979) used the extraction term proposed by Gardner and Ehlig
(equation 2.20) in a mechanistic simulation model describing the process of soil moisture
extraction by root systems. They conducted experiments to test the model's validity and
found that root resistance was the dominant resistance term through most of the soil moisture
availability range. The increase in soil hydraulic resistance resulting from soil moisture
extraction did not become significant relative to root resistance until the soil moisture reserve
was substantially depleted.
Nimah and Hanks (1973a) developed a model for estimating soil water, plant, and
atmospheric interrelations. Their onedimensional model predicted water content profiles,
root extraction, and root water potential under field conditions. It used a root extraction
term of the form
[Hot + (1.05z) h(z,t) s(z,t)] RDF(z)K(0)
S(z,t) = (2.22)
AxAz
where
Hroot = the effective water potential in the root at the soil surface,
h = soil matric potential,
s = osmotic potential,
RDF = proportion of total active roots in depth increment Az,
K hydraulic conductivity, and
Ax = distance between the plant roots at the point where h and s are measured.
The model predicted significant changes in root extraction due to variations in
pressure headwater content relations and root depth. A field test of the model (Nimah and
Hanks, 1973b) with an alfalfa crop provided best agreement between predicted and computed
values 48 hours after water addition. Poorest agreement was found immediately following
irrigation.
Feddes et al. (1976) defined the extraction term in two different ways. The first
extraction term was dependent on the hydraulic conductivity of the soil, the difference in
pressure head between the soil and the rootsoil interface, and a root effectiveness function.
The second term was defined as a function of soil water content. The governing partial
differential equations were solved by a finite difference and a finite element approach, and
the computed values were verified against field measurements. Only slight discrepancies
between the two numerical methods were found and these were attributed to the different
manner in which they handle boundary conditions and water uptake by roots. Agreement
between measured and computed values ranged from excellent to satisfactory.
Slack et al. (1977) developed a mathematical model of water extraction by plant roots
as a function of leaf and soil water movement and water uptake. A twodimensional radial
flow model was used to describe water movement and estimate transpiration from corn grown
in a controlled environment under soil drying conditions. A microscopic and a macroscopic
approach was used to produce a combined extraction model which predicted daily
transpiration adequately for the period modeled.
Feddes et al. (1978) developed an implicit finite difference model to describe water
flow and extraction in a nonhomogenous profile under the influence of groundwater. Their
extraction term was considered to be a function of the soil moisture pressure head, ti. The
extraction rate was held at a maximum level if was above a set limit, 0. When the pressure
head fell below 0i, water uptake decreased linearly with / until the wilting point, #3, was
reached. Below 0i water uptake by roots was set at zero. The extraction term between iz and
0i was defined by
S() = Sma (2.23)
S 3
where
Sma = the maximum possible transpiration rate divided by the effective rooting
depth.
The model yielded satisfactory results in predicting both cumulative transpiration and
distribution of soil moisture content with depth.
Tollner and Molz (1983) developed a physically based extraction function for use in
the general flow equation. The function assumes that water uptake per unit volume of soil
is proportional to the product of contact length per unit soil volume, root permeability per
unit length, and water potential difference between soil and root xylem potential. A factor
which accounts for reduced rootsoilwater contact as water is removed was included in the
extraction term. Model predictions compared favorably with experimental results from the
literature.
Studies conducted by Denmead and Shaw (1962) and Ritchie (1973) with field crops
showed that during periods of high potential evapotranspiration (ET) demand, actual ET was
considerably less than the potential rate even though the available soil water was considered
adequate. They also observed that during periods of low potential ET demand actual ET
equaled potential rates, even down to very low soil water contents.
Smajstrla (1982) developed an extraction function that was consistent with these
observations. His function was based on potential extraction rate and soil water status and
was defined as
S = TpR (T/CR) (2.24)
where
Tp = potential transpiration rate, and
C = a calibration constant.
R was defined as the relative available soil water content and was given mathematically by
R (2.25)
8t. 8 p
where
S= soil water content,
0P, = soil water content at the wilting point, and
re = soil water content at field capacity.
The model was calibrated using the rate of soil water content change observed during
lysimeter field experiments with soybeans. The function permits a rapid rate of soil water
extraction when soil water is readily available. As the water is depleted, the extraction rate
declines logarithmically. Extraction is limited during periods of high potential ET and is
allowed to recover to near potential rates during periods of low potential ET. The function
limits extraction to very low rates as the permanent wilting point is approached.
Stone (1987) modified Smajstrla's extraction term to include a root distribution
function and substituted potential transpiration (Tp) with actual ET, thus eliminating the
calibration constant. The extraction term was used with success in a twodimensional model
simulating extraction patterns under agricultural crops.
Lvsimetrv
Lysimeters have been used for determining soilplantwater relationships for many
years. A review of lysimetry was published by Aboukhaled et al. (1982). Lysimeters can
generally be described as weighing and nonweighing (or drainage) types whose advantages
and disadvantages were discussed by Smajstria (1985). He found that weighing lysimeters
are more costly and difficult to maintain because of requirements for service and calibration.
The tanks of weighing lysimeters must also be constructed with sufficient mechanical strength
to support the soil profile without deformation and be supported by expensive load cells.
They are ideally suited to studies in which hourly or daily soil water budgets are required.
Drainage lysimeters are cheaper to construct because the surrounding soil can be used
to support the tanks (Aboukhaled et al., 1982). They are used for studies in which weekly
or monthly soil water budgets are of interest because instrumentation to measure soil water
status is not sufficiently accurate. Smajstrla (1985) developed a set of drainagetype field
lysimeters for water management studies on deep sandy soils for Florida humid climatic
conditions.
Feddes et al. (1978) studied the influence of the groundwater table on the amount of
water available for transpiration of red cabbage on sticky clay with a drainagetype field
lysimeter. Their system was designed to continuously maintain the water table in the
lysimeter at approximately the same depth as it was in the surrounding field. Drainagetype
lysimeters with water table controls for Florida conditions were described by
Clayton et al. (1942), Speir (1962), and Shih et al. (1977).
CHAPTER 3
EXPERIMENTAL PROCEDURE
The Field Lysimeter System
A field lysimeter system was constructed to study soil water redistribution and
extraction patterns of three different soil profiles irrigated by a line source under the
influence of a shallow water table. The system consisted of 6 rainsheltered drainagetype
lysimeters arranged in tandem to simulate rowcrop conditions. Rain shelters were used to
eliminate unwanted rainfall from the lysimeters during water management studies
(Figure 3.1). Soil water status data were collected with a microcomputerbased data
acquisition system which consisted of tensiometermounted pressure transducers, a data
acquisition and control device, and a portable microcomputer.
Lvsimeter Construction
Each of the six lysimeters constructed for this facility provided a rectangular cross
section of 0.7 m2 and an effective surface production area of I m2. Their profiles were
0.75 m deep.
The lysimeters were built from commercially available galvanized steel
3 x 0.9 x 0.75 m water troughs. Each trough was partitioned to create two lysimeters. The
partitions were machined from sheets of aluminum and bonded into the troughs with an
epoxyresinbased metal filler. The strength and water tightness of each bond was tested by
filling one of the lysimeters in each trough with water while leaving the other empty and
inspecting for leaks. The aluminum was then covered with a primer, an epoxybased paint,
and an asphaltbased paint to prevent corrosion.
The lysimeters were installed at the Institute of Food and Agricultural Sciences Irrigation
Research and Education Park at Gainesville, Florida. A trench was excavated to a depth of
25
z0
LUHI
FF
00U)
0 UL
ii 4)
I, C) W
z~J 0
~, m
LL=
LUi IJi Z 0 IC
LrC j II
D>II C
m0 0 w <)
i t LO l0
Ul 31: IM
II <
2 I 0
00
ma
U) m 0 0 0 '1
0 0
ogD I L ) 0
oo CI I
III
I 2M
Fr L U
O LU I
LiL
/ 0
0 0t
0/ 1/7 L I>
llz
LLCl ;
iLwu Li
M )L
m n 0
cf)
0.5 m with a backhoe and the remaining 0.25 m were excavated by hand to ensure that the
bottom surface of the trench was level. A thin layer of coarse builder's sand was placed in
the bottom of each trench, leveled, and tamped to ensure that the bottoms of the lysimeters
were completely supported. The lysimeters were installed in tandem, with a 0.3 m spacing
between troughs. A carpenter's level was used to ensure that they were level in order to
avoid nonuniform water table heights within the lysimeters. A standpipe and a drainage/
distribution system were placed in each of the lysimeters, and a 50 mm layer of builder's
sand was laid over the distribution systems to facilitate drainage and water table control. Soil
was backfilled around the lysimeters while soil profiles were simultaneously reconstructed
in the lysimeters by hand packing to prevent the shapes of the lysimeters from distorting
while they were being filled.
Three different soils were used to fill the lysimeters, with each of the soils occupying
a pair of lysimeters. Because Florida agricultural soils are singlegrained fine sands, the soils
used in the lysimeters were chosen accordingly and further selected to provide high,
intermediate, and low waterholding capacities. The soils used were Pomona sand, Astatula
fine sand, and Builder's sand, respectively.
Profile reconstruction. Handpacking was used to assure that the individual soil
profiles were as homogeneous as possible (Smajstrla, 1985). Lysimeter profiles were
reconstructed in approximately 50 mm thick layers by a continual process of filling, tamping,
and raking between layers.
To simulate row crop conditions, beds were formed within the lysimeters. One 0.9 m
long trapezoidal bed ran across the width of each lysimeter as illustrated in Figures 3.1 and
3.2. The soil profiles were reconstructed uniformly for the first 450 mm and then the beds
were shaped. Each bed was 1.0 m wide at its base, 300 mm wide on top, 300 mm tall, and
side slope ratios were 3:2. Half the height of the bed (150 mm) protruded above the lip of
the lysimeter. To contain the ends of the beds, sheets of galvanized steel, cut into trapezoidal
shape, were pushed into the soil profile along the inner edges of the lysimeters. The
irrigation laterals were then placed and the beds covered with black plastic mulch.
Lvsimeter water table control. A system was designed to provide accurate regulation
of the water table within the lysimeters. The design consisted of a float valve and a sump
which were located in a standpipe.
The standpipes were constructed by cutting a 250 mm diameter polyvinyl chloride
(PVC) pipe into lengths of 650 mm. Holes (3 mm) were drilled into the bottom 100 mm of
each pipe to facilitate the flow of water. Very fine mesh nylon screen was used to cover the
holes and prevent soil particle loss. The standpipes were then glued to the bottom of the
lysimeters in an upright position (Figure 3.2). Lids were manufactured to keep out debris and
insects.
Attached to the standpipe was a distribution system (Figures 3.1 and 3.2) to facilitate
the movement of water along the bottom of the lysimeter. PVC pipe elbows were used to
connect lengths of 26 mm well screen to form a rectangular distribution system. One corner
of the distribution system led into the standpipe. On the opposite end of the distribution
system, a reducing elbow was used to install a vertical 16 mm PVC pipe and a gas jet which
allowed water to be extracted to the surface by vacuum. The gas jet was used for rapidly
extracting water or for lowering the water table below the level allowed by the sump.
The water table in each lysimeter was recharged from the float valve (Figure 3.3)
which consisted of a water feeder valve regulated by a float attached to the valve lever arm.
A coupling and slip sleeve were used to connect 6.4 mm diameter copper tubing to the float
valve. Tygon tubing was slipped over the copper tubing and led to the top of the standpipe.
From there, 10 mm polyethylene (PE) tubing led to a 40 L water supply tank which provided
the positive head required by the float valve. Two nuts, welded together, were used to mount
the float valve onto a threaded stand and permitted fine adjustment of the water levels
(Figure 3.4). A minimum level of 20 mm was required for the float valve to operate properly.
Once calibrated, the float valve allowed the water level to be adjusted to within 3 mm of the
desired depth setting. The volume of water added to each water table through a float valve
was determined by measuring the volume of water depleted from the corresponding water
tank. To accurately measure depletions, the water table tanks were replenished with the use
of a graduated cylinder.
Figure 3.3. Float valve and sump used to control the water table level.
Figure 3.4. End view of the float valve showing details of the mechanism
used to mount the valve onto a threaded stand.
The sumps were constructed to eliminate fluctuations in the water table caused by
percolating irrigation water. Any water draining into the sumps was extracted by vacuum
to individual drain traps. The sumps were made by gluing a bottom to one end of a 30 mm
length of 50 mm PVC pipe. Each sump was then attached to the lysimeter bottom within the
standpipe (Figures 3.2 and 3.3). Vacuum tubing, which led to individual drain traps, was
connected to the sumps with a hose barb.
A portable drain trap was constructed for each lysimeter from 100 mm PVC pipe.
PVC endcups were then glued on a 600 mm length of pipe to form a cylindrical tank with
a capacity of approximately 5 L. A graduated cylinder was used to accurately measure the
volumes of water collected in the drain trap.
The vacuum system consisted of a 250 L/min vacuum pump, a 300 L vacuum tank,
a regulator to start and stop the vacuum pump, and vacuum control valves to permit the
desired vacuum to be set.
Irrigation system. The irrigation system was a low flow, manually operated line
source drip irrigation system. Each lysimeter was irrigated by a separate lateral, which was
supplied from a 60 L water tank. The tank, placed on a platform, provided the positive head
required by the irrigation system. Laterals were placed along the centerlines of each bed
(Figures. 3.1 and 3.2). Twelve drip emitters, closely spaced at 75 mm, in order to simulate
a line source, were installed on each lateral. Spot Vortex emitters were used with flow rates
of 2.4 L/hr or 32 L/hr per m lateral length. To accurately measure the irrigation volumes,
a graduated cylinder was used to add water to the irrigation tanks.
Rain shelter. An automated rain shelter was used to cover the lysimeters in the event
of rainfall. The system consisted of a shelter mounted on a rail system, electric motors and
controls, and a rainfall detector and switch to activate the motors. The shelter was part of
a preexisting lysimeter system (Smajstria, 1985).
The Data Acquisition System
A microcomputerbased data acquisition system for continuous soil water potential
measurements was developed. The system consisted of tensiometermounted pressure
transducers, a Keithley System 570 work station data acquisition and control device, and a
Zenith 161 portable microcomputer.
Pressure transducers were mounted on Model S Irrometer tensiometers by drilling
3 mm holes into the tensiometer reservoirs opposite the vacuum gage ports. The pressure
transducers used were Micro Switch model 141PC15D transducers which measure pressures
from 0 to 100 kPa. Each transducer produced an analog output of I to 6 volts which was
proportional to the vacuum exerted on a membrane in the transducer. One side of the
membrane was in contact with the tensiometer water while the other side was open to the
atmosphere. As pressure on the water changed, the membrane deformed proportionally, and
the voltage across the membrane changed. The pressure transducers required an 8volt direct
current (VDC) power supply at 8 mA each. These transducers were used because they were
shown to be linear, temperature compensating, and adaptable to field use (Stone et al., 1985).
The System 570 data acquisition system had provisions for a maximum of 32 single
ended input channels. To provide flexibility, 12 conductor signal cable was installed, thereby
providing 12 available signal leads for each pair of lysimeters. The signal cable was 22
American Wire Gage (AWG) with 6 pairs of stranded chrome conductors, each pair shielded
to prevent induced voltage errors. A 12gage power supply cable was installed to provide
power to a junction box for each pair of lysimeters.
The microcomputer and the data acquisition system were housed in an air conditioned
instrumentation building about 60 m from the lysimeters. Between the instrumentation
building and the junction boxes, 3 signal cables and 3 power supply cables were buried in
a 75 mm PVC pipe. One data cable and one power supply cable led into each junction box
and a three conductor cable connected each transducer to its junction box. Two of the
conductors were used for power supply, while the third carried the signal from the
transducer.
Three of the 6 lysimeters were instrumented with 10 tensiometers per lysimeter,
engaging 30 channels (Figure 3.1). One channel was used to monitor atmospheric pressure
and the remaining channel was used as a backup. In order to provide a complete description
of the soil profile in the lysimeter, the 10 tensiometers in each of the instrumented lysimeters
were installed at depths and distributions as follows: 3 tensiometers were installed at each
of 3 depths, 100 mm, 200 mm, 300 mm, and at distances of 50 mm, 150 mm, and 250 mm
from the lateral. One tensiometer was installed at a depth of 500 mm (Figures 3.5 and 3.6).
The Keithley System 570 data acquisition and control device was the interface
between the microcomputer and the tensiometer pressure transducers. It consisted of a
chassis/mother board and an interface card for the microcomputer. The mother board had
individually numbered screw terminals to which the pressure transducer data cables were
attached.
The System 570 was supported by the Soft500 software package which is an extension
of Advanced BASIC and provided foreground/background architecture, multitasking, array
management, memory management, disk access and storage, and a library of input and output
commands accessible directly from BASIC. Real time analog measurement was accomplished
by writing a Basic program incorporating Soft500 commands.
Data Acquisition System Operation. The data acquisition system permitted soil matric
potentials to be monitored nearly continuously. Data from the 30 tensiometerpressure
transducers were collected at one minute intervals and averaged over a 10minute period.
The average soil matric potential of each tensiometer for that period and the corresponding
time were saved to disk. The 320 kbyte 5 1/4 inch floppy disks used for storage of the soil
water potential data permitted 7 days of continuous data collection.
Instrument calibration. The 30 pressure transducers were individually calibrated to
analyze their response characteristics. Transducer outputs were measured at decreasing
pressures from 0 to 40 kPa and individual transducer calibration curves were produced.
Regression analysis of the data indicated that the calibration curves were linear with R2
values of 1.000.
RADIAL DISTANCE FROM DRIP LATERAL (mm) 3
250 150 50 0 150 300
/ 200 \
I/ \
...._ g g 300 .
500
E
Figure 3.5. Side view of a lysimeter showing the location of the
10 tensiometers.
100
200
S300
E
S500

SPACING ALONG BED Cmm3
150 450
750 900
f5u
Figure 3.6. End view of a lysimeter showing the location of the
10 tensiometers.
 
*
I
Calibration curves of the assembled tensiometerpressure transducer systems were also
developed to take into account the effect of the hanging tensiometer column on the
transducer output. The curves were developed by exerting a known vacuum on the
tensiometer ceramic cups and reading the transducer outputs. The curves again proved to
be linear with R2 values ranging from 1.000 to 0.998. The net effect of the hanging column
was to shift the output voltage of the transducer upwards by a an amount directly
proportional to the column length. Figure 3.7 shows two typical calibration curves of
transducer output voltage versus applied pressure. The lower curve shows the data and
regression line of an individual transducer. The upper curve is the output of the same
transducer mounted on an 18 inch (457 mm) tensiometer.
Pr terssre ssob)
4
acquisition system provided 12bit resolution. This produced 4096 counts in the 05 volt
5 10 15 ?0 25 30 35 40
Pressure kPo)
Figure 3.7. Calibration curves of output voltage versus negative
pressure for transducer no. 4.
Pressure measurement accuracy. The analogtodigital (A/D) board in the data
acquisition system provided 12bit resolution. This produced 4096 counts in the 05 volt
range, 0.0012 volts per count, or a resolution of approximately 0.03 kPa. The relationship
between voltage and pressure for transducer no. 4 is given by equation (3.1).
kPa = (volts 1.1153) / 0.04766
35
The above equation resulted in a maximum expected error of 0.06 kPa at 1 volt
(0.3 kPa of pressure) and 0.27 kPa at 5 volts (83.7 kPa of pressure) (Stone et al., 1985).
Stray voltages in the instrumentation building reduced resolution to 0.01 volts or 0.2 kPa.
Instrumentation error of 0.2 kPa produces error of 20% to 4% in pressures from 1 kPa to
5 kPa, respectively. Therefore, it was important to recognize that significant errors could
result when measuring soil water potentials near saturation.
Atmospheric pressure effects. Soil moisture matric potentials were measured using
differential pressure transducers with atmospheric pressure as their reference. To determine
the effects of atmospheric pressure changes on the performance of the differential pressure
transducers, ambient pressure was measured using an absolute pressure transducer. It was
found that atmospheric pressure did not vary enough to warrant correction of the matric
potential measurements. Typical fluctuations varied from 0.1 to 0.3 kPa. These were well
within the range of other variations in measurement accuracy. Figure 3.8 illustrates typical
variation of ambient air pressure during a week in June, 1988.
0.50
0j .25
        
a0.
JUNE
Daily S
13 14 15 16 17 18 19
JUNE
Daily Scale
Figure 3.8. Fluctuations in atmospheric pressure for the week of
June 1319.
Lvsimeter System Operation
In Florida, vegetables are typically grown under conditions similar to those that were
simulated in the lysimeters. The tomato was chosen for this study for several reasons:
1. Tomatoes are one of the most important cash crops in the state.
2. Tomato plants have high transpiration rates which would allow measurable
extraction rates and produce very definite redistribution and extraction patterns.
3. Information on redistribution and extraction patterns, contribution of the water
table to crop water use, and water budgets was critically needed.
The Super Bush (VFN) hybrid (Lycopersicon esculentum Mill.) was selected for its
vigorous growth and large, leafy bush. These were desirable traits because they increase the
transpiration rate of the plant. Seedlings approximately 100 mm tall were purchased from
a commercial nursery and transplanted to the lysimeters on March 15, 1988. Two seedlings
were planted along the centerline of each bed, 450 mm apart, and each 225 mm from the
edge of the lysimeter. Plant growth was monitored by taking biweekly canopy height and
width measurements with a meter stick. Mature plants were removed from the lysimeters on
June 26, 1988 which allowed for a three month growing season.
Irrigation scheduling. Three of the six lysimeters were instrumented. The two
leftmost lysimeters in Figures 3.1 and 3.9 contained the Astatula f.s. and were labeled El and
E2 (for East I and East 2). The two center lysimeters contained the Builder's sand and were
labeled Cl and C2. Finally, the two rightmost lysimeters contained the Pomona sand and
were labeled Wl and W2. The three instrumented lysimeters were E2, CI, and WI.
Astatula f.s. Builder's sand Pomona sand
El IE C3 C2 ul W E
Figure 3.9. Arrangement of the three soils in the lysimeters.
An irrigation was applied when the soil water potential in a lysimeter dropped below
10 kPa. Soil water potential was used as the trigger point for irrigation instead of a
percentage of water depletion from field capacity because potential is independent of soil
type. Thus, the plants were subjected to a uniform stress level regardless of soil type. The
10 kPa value was chosen because it is consistent with the current practices used by growers
who tend to keep the root zone moist and the plants stressfree.
The irrigation system was manually operated and the lysimeter system was serviced
daily. Every morning the transducer outputs were monitored. If any of the instruments read
below 10 kPa, the appropriate lysimeters were irrigated.
Permanent ink was used to mark reference levels on the irrigation and water table
tanks. The reference level in the irrigation tanks was placed at the 40 L mark, or at 3/4 full.
To provide the positive head required by the irrigation system, the irrigation tanks were
always maintained at the reference level. The irrigation process began by adding the desired
irrigation volume, carefully measured with a graduated cylinder, to the water already in the
irrigation tanks. The valve on the irrigation tank was then opened and the water level in the
tank allowed to drop to the reference level.
Immediately following irrigation with the drip emitters, the volume of water used to
maintain the water table at the set level since the previous irrigation event was determined.
Initially, the water table tanks were filled to the reference level (40 L mark). As the water
table was recharged, the water level in the tanks dropped. After each irrigation event, water
was added to the corresponding water table tank with a graduated cylinder until the tank was
refilled to the reference level.
The vacuum system was left running continuously and the drain traps were checked
daily. If any water had collected in the drain traps, it was measured and recorded.
Irrigation volumes. The objective in determining the irrigation volumes to be applied
was that they should be large enough to wet the soil profile to field capacity but not so large
that they would drain from the root zone. The approximate irrigation volumes for each soil
type were estimated using the soilmoisture retention curves. The actual volumes were
determined through trial and error by incrementing the estimated volumes until drainage just
began. Adjustments were made to the volumes over the growing season, because as the plants
grew, they extracted larger volumes of soil water which had to be replenished.
Table 3.1 presents the irrigation volumes for the mature plants and their
corresponding depths over the surface area of the lysimeters. These data demonstrated the
need to apply small, frequent irrigations to minimize drainage. They also demonstrated that
the volumes were a function of soil type and that the soils with low waterholding capacities
(the Builder's sand and the Astatula f.s.) required much smaller irrigations than the Pomona
sand. However, soil type does not necessarily influence the plant water requirement.
Table 3.1. Volumes of Irrigation Applications and
Equivalent Depths Over the Lysimeter Areas.
Soil type Volumes (L) Depths (mm)
Builder's sand 4.00 3.46
Astatula f.s. 4.00 3.46
Pomona sand 8.00 6.93
Irrigation began immediately after the seedlings were transplanted but the volumes
were not recorded until the data acquisition system became operational on April 1. By
April 15, a water table 720 mm from the surface had been established in all the lysimeters
and was being maintained by the float valves. On this date, collection of data on the
contribution of the water table to crop water use began.
Soil Characteristics
Samples of each soil were taken and submitted for analysis to the University of
Florida's IFAS Soil Testing Laboratory. Analyses were done on pH, percent organic matter,
macro and micro nutrient concentrations, copper toxicity, and soluble salts. Organic matter
and pH data are included in Table 3.2. Complete results can be found in Appendix A.
39
Analysis of the Pomona sand indicated that the soil had an initial pH of 4.6. Because
the recommended soil pH for tomatoes is 6.06.8 (Stephens, 1984), 7.8 Mg/ha of dolomitic
limestone was incorporated into the surface 25 cm of soil to adjust the pH level, as
recommended by the Soil Testing Laboratory. The liming recommendation was based on the
AdamsEvans lime requirement test.
Table 3.2. Characteristic Soil Properties.
Bulk Organic Saturated Residual Saturated
Soil density matter water cont. water cont. Hydraulic Cond. pH
type (g/cm ) (g/g) (mm3/mm3) (mms/mms) (mm/hr)
Pb % r0 0r K,
Builder's sand 1.68 0.01 0.35 0.01 536.07 6.8
Astatula f.s. 1.66 0.58 0.41 0.03 428.35 7.1
Pomona sand 1.10 5.34 0.54 0.08 169.77 4.6
The soil test reports also provided standard fertilization recommendations for each
soil which were used to compute the appropriate fertilizer applications. The
recommendations were similar for all three soils and called for 135 kg N, 180 kg P2Os, and
180 kg K20 per ha. Half the requirement was met by mixing a complete fertilizer with a
micronutrient supplement into the soil prior to planting. The remaining requirement was met
by adding a complete water soluble fertilizer to the irrigation water tanks and applying it
through the irrigation system at each irrigation. This fertilizer was supplemented with a
micronutrient foliar spray.
Bulk Density
The bulk densities of the three soils used in the lysimeters were determined in the
laboratory by packing five samples of each soil into fifteen aluminum cans which were
50 mm deep and had a 50 mm diameter. After oven drying at 105" C, the mass of each
sample was determined by subtracting the mass of the can from the total mass. The can
volumes were found by determining the mass of water required to completely fill each can
and converting that mass to volume. Excellent agreement was obtained between the five
samples of each soil which were averaged to produce the bulk densities given in Table 3.2.
SoilMoisture Retention Curves
Soilmoisture retention curves were developed for the three soils used in the study
(Figure 3.10). A pressure plate extractor and the previously determined bulk densities were
used to develop desorption curves of volumetric soil water content versus pressure potential.
Five soil samples were used to determine the volumetric water content of the soil at a given
pressure. Early results indicated that the time required for soil samples to reach equilibrium
at a given pressure was much larger than the 24 hrs recommended in the literature. In this
study, a five day period was used to ensure equilibrium. Saturated and residual water
contents are given in Table 3.2.
The coefficients of variation of the water contents near saturation (5 to 1 kPa)
ranged from 2% to 8%, respectively, emphasizing the difficulty in obtaining these data in this
pressure range. Some of this variability may have been due to normal sampling variability,
while much was probably due to variation in control pressure in this range. The coefficients
of variation generally decreased to less than 0.1% as potential decreased to 100 kPa.
Saturated Hydraulic Conductivity
Saturated hydraulic conductivity (K,) was determined with a constant head
permeameter. A soil core of length AL was packed into a permeameter and saturated with
water from the bottom upwards. A constant depth of water was then established over the soil
column and the system allowed to reach equilibrium. The outflow of water from the system,
Q, was measured at fixed time intervals, t, to determine when equilibrium conditions were
established. The constant head difference, At, was calculated by adding the depth of the
standing water to AL. The saturated hydraulic conductivity was determined from
Q AL
K, A (3.2)
At AO
where A is the crosssectional area of the soil column. The saturated hydraulic conductivities
of the three soils are given in Table 3.2.
Hydraulic Conductivity Curve
Hydraulic conductivity curves were developed for each of the three soils by the van
Genuchten (1980) method because it has been extensively used to successfully predict
conductivity curves of sandy soils. The observed data required by the method for each soil
are the soilmoisture retention curve, the saturated and residual water contents, and the
saturated hydraulic conductivity.
First, equation (2.6) was fitted to the data. The resulting parameters a, n, and m are
tabulated below. To test the fit of the analytical functions, the predicted soilmoisture
retention curves were plotted against the observed curves in Figure 3.11. The fits were
excellent for the builder's sand and the Astatula f.s. and good for the Pomona sand.
Equation (2.9) was then used to predict hydraulic conductivity which was plotted against
pressure head in Figure 3.12.
Table 3.3. Parameters for van Genuchten's Equation.
Soil type a n m
Builder's sand 0.032 6.59 0.848
Astatula f.s. 0.033 3.53 0.717
Pomona sand 0.029 3.01 0.668
O Builder's sand
T Astatula f.s.
E Pomona sand
Si
7 1000
Q)
U)
1000
1 0 0 .. . . . . .
0.00 0.10 0.20 0.30 0.40 0.50
Soil Moisture Content, 6
Figure 3.10. Soilmoisture retention curves for the three soils used in
this study.
43
1E+04
Observed curves
0  Fitted curves
SCurves from left to right:
E Builder's sand
E Astatula f.s.
S \ Pomona sand
U
(D
T 1000
J)
0.00 0.10 0.20 0.30 0.40 0.50
Soil Moisture Content, 0
Figure 3.11. Fitted and observed soilmoisture retention curves for
the three soils used in this study.
102  .uu UII 0 oZ UIIU _
S  Astatulo f.s.
\  Pomona sand
E \ \ \
E
10 2
U
._ 106 \\
0 104
0
\ 1
10o
S10 1I ,
10 10
102 103 1 0
Negative Pressure Head (mm)
Figure 3.12. Hydraulic conductivity curves for the three soils used in
this study.
CHAPTER 4
EXPERIMENTAL RESULTS
Plant Growth
The months of March and April, 1988, were unusually cool in Gainesville.
Consequently, plant growth and water use were uncharacteristicly limited during those
months. The majority of plant growth occurred in May and June (Table 4.1). For this
reason, the discussion of water use was limited to these final eight weeks of the season.
Table 4.1. Canopy Dimensions (m).
Soil Type March 15 May 01 June 26
and
Lysimeter Height Width Height Width Height Width
Builder's Sand Cl 0.11 0.06 0.24 0.30 0.70 0.45
Builder's Sand C2 0.10 0.07 0.26 0.28 0.67 0.38
Astatula f.s. El 0.12 0.05 0.25 0.31 0.73 0.55
Astatula f.s. E2 0.10 0.07 0.27 0.26 0.86 0.48
Pomona sand W1 0.10 0.07 0.26 0.29 0.93 1.06
Pomona sand W2 0.11 0.06 0.28 0.32 0.93 1.14
Despite efforts to maintain plant size uniform across all soil types, the plants in the
Pomona sand grew considerably better than those in the Astatula f.s., while the smallest
plants were produced in the builder's sand. Yield data were not collected because the fruit
were afflicted by blossom end rot, particularly those in the Pomona sand. The disease, which
is caused by calcium deficiency, was treated by spraying the plants with a foliar calcium
supplement. The plants responded well to the treatment, but not quickly enough to produce
a representative yield by the end of the experiment.
Measurement Accuracy
There are inherent limitations when drainage type lysimeters are used to conduct
water use studies. These limitations which, along with other measurement errors, may reduce
the accuracy of collected data, must be recognized and accounted for.
Water Budget Accuracy
All water volumes applied to or extracted from the lysimeters were measured with
a graduated cylinder. The accuracy of the graduated cylinder and the ability to refill the
water supply tanks reproducibly determined the accuracy of the measured water volumes.
A I L graduated cylinder with accuracy to 10 mL was used to refill the supply tanks and
measure drainage water extractions. Calibration tests indicated that the supply tanks were
refilled reproducibly to within 50 mL which produced errors in water inflows and outflows
of less than 2%.
Weekly. Monthly, and Seasonal Water Use Accuracy
For accurate shortterm measurement of water use from a lysimeter, an accurate
measurement of soilwater storage or soilmoisture content is required. Because drainage
type lysimeters do not permit a direct determination of soilmoisture content, it was
determined indirectly by monitoring soil water potential.
Soilmoisture content and water potential are related through the soilmoisture
retention curve which was measured for each soil. But, as discussed earlier, the translation
of soil water potential to soilmoisture content on the basis of the soilmoisture retention
curves measured introduced errors up to 8% near saturation. The combination of
instrumentation error of 0.2 kPa (20% to 4% of I to 5 kPa) and error resulting from the
pressure extractor procedure were the causes of this error. The resulting error prevented
accurate measurements of daily water use, although wetting or drying patterns were still
determined. As the measurement period increased, the relative size of the error diminished
with respect to the water used by the crop. Thus, weekly water use was accurately
determined, especially during high water use periods.
Crop Water Use
Because the research was conducted in rainsheltered lysimeters buffered by well
watered grass, the water use data may not represent the water use requirements of an
extensive tomato crop grown under field conditions in Gainesville, Florida. The data are
primarily a comparison of crop water use contributions by the water table and the irrigation
system with respect to soil type.
Partitioning of Crop Water Use
For the purposes of this study, measured crop water use was defined as
crop water use = irrigation drainage + water table recharge (4.1)
The design of the lysimeters permitted crop water use to be partitioned into irrigation and
water table recharge components. Monthly values in terms of mm over the surface area of
the lysimeter for each of the components are presented in Figures 4.1 and 4.2. The numbers
over the irrigation and water table recharge bars represent the percentage of crop water use
contributed by each component. The number over the drainage bar represents the percentage
of irrigation water lost to drainage. Drainage was a function of soil type. The Pomona sand,
with the greatest water holding capacity, consistently had the smallest drainage fraction.
Overall, drainage volumes were small because, as discussed earlier, irrigation applications
were designed to minimize deep percolation.
During May, when the plants were still maturing, the contribution of the water table
to crop water use in the builder's sand, the Astatula f.s., and the Pomona sand was 16%, 19%,
and 33%, respectively (Figure 4.1). In June, when the plants were mature, the water table
contributions were 26%, 34%, and 31% (Figure 4.2). These data indicate that as the plant
root systems developed and extended down towards the capillary fringe in the builder's sand
and the Astatula f.s., the water table contribution increased significantly. In the Pomona
sand, on the other hand, plant size had little effect on water table contribution. In fact, the
contribution of the water table to crop water use in Pomona sand remained the same
PARTITIONING OF CROP WATER USE
May, 1988
Depth. (mm)
Buider's and Asttull f I. Pomona sand
E Crop water use Irom Irrigation M Crop wtllr ul fhror watr ftble
E3 Orainago
Figure 4.1. Partitioning of crop water use in May.
PARTITIONING OF CROP WATER USE
June, 1988
Deplhs (mm)
142%
Builder'l land Astatula I.e. Pormon sand
E[ Crop waler usl from Irrigation S Crop water usl from water table
q Drainage
Figure 4.2. Partitioning of crop water use in June.
throughout the growing season. Consequently, it was concluded that the properties of the
soil, and in particular the particle size distribution and organic matter content, allowed the
capillary fringe to rise high enough in the Pomona sand to make water available to the plant
roots, regardless of plant size.
It is also possible that the rate of growth of the plant root system towards the capillary
fringe in the builder's sand and the Astatula f.s. was a function of the irrigation schedule.
If these two soils could not hold the daily crop water requirement from only one irrigation
event per day, the root system may have developed towards the water table in response to the
water available there. Pressure transducer readings of up to 50 kPa indicate that the plants
in all three soils may have experienced mild, daily short term water stresses during the hours
preceding irrigation. At no time, however, did the plants show any visual signs of water
stress such as wilting.
Weekly, Monthly, and Seasonal Water Use
Weekly crop water use for the final eight weeks of the growing season is presented
in Figures 4.34.8. Each figure contains two graphs. The graph on the left shows weekly
water use in the form of stacked bars where the lower bar depicts the contribution of
irrigation water to crop water use and the upper bar the contribution of the water table. The
other half of each figure presents potential evapotranspiration (computed by the Penman
equation) over the same period. By comparing the two graphs, it can be clearly seen that the
fluctuations in water use correspond closely with the fluctuations in potential ET.
The water use data indicate that crop water use was a function of plant size and that
the plants' source of water (irrigation versus water table) was a function of soil type. The
plant size data presented in Table 4.1 do not adequately illustrate the striking difference in
appearance between the plants. This is more apparent in Figures 4.14.8 where it can be
clearly seen that although the water use patterns were the same for all three soils, the total
water use values varied dramatically from soil to soil.
WEEKLY WATER USE
Depth (mm) MAY
60
50
40
30
20
10
0;
*B
08 15 22
Weekly Scale
From irrigation
From water table
WEEKLY POTENTIAL ET
Depth (mm) MAY
60
50
40
30
20
10
08 15 22 2S
Weekly Scale
Potential ET
Figure 4.3. Weekly crop water use in builder's sand, May, 1988.
WEEKLY WATER USE
Depth (mm) JU
60
50
40
30
20
10
]a
NE
05 12 19
Weekly Scale
From irrigation
From water table
WEEKLY POTENTIAL ET
Depth (mm) JUNE
60
50
40
30
20
10
0
05 12 19 21
Weekly Scale
Potential ET
Figure 4.4. Weekly crop water use in builder's sand, June, 1988.
WEEKLY WATER USE
Depth (mm) MAY
60
50
40
30
20
10
08 15 22 29
Weekly Scale
From irrigation
] From water table
WEEKLY POTENTIAL ET
Depth (mm) MAY
60
50 .   
0
30
20
10
0
08 15 22 2!
Weekly Scale
Potential ET
Figure 4.5. Weekly crop water use in Astatula f.s. May, 1988.
WEEKLY WATER USE
Depth (mm) J
JUNE
60
50
40
30
20
10
05 12 19 26
Weekly Scale
From irrigation
f From water table
WEEKLY POTENTIAL ET
Depth (mm)
JUNE
60
50
40 " ''
30
20
10
0
05 12 19 21
Weekly Scale
Potential ET
Figure 4.6. Weekly crop water use in Astatula f.s. June, 1988.
WEEKLY WATER USE
Depth (mm) MAY
60
50
40
30
20
10
0
08 15 22 29
Weekly Scale
E From irrigation
E From water table
WEEKLY POTENTIAL ET
Depth (mm) MAY
60
50
40
30
20
10
0
08 15 22 29
Weekly Scale
Potential ET
Figure 4.7. Weekly crop water use in Pomona sand, May, 1988.
WEEKLY WATER USE
Depth (mm)
JUNE
60
50
40
30
20
10
0 
05 12 19 26
Weekly Scale
From irrigation
[] From water table
WEEKLY POTENTIAL ET
Depth (mm)JUNE
JUNE
60
50
40
30
20
10
05 12 19 2E
Weekly Scale
Potential ET
Figure 4.8. Weekly crop water use in Pomona sand, June, 1988.
Seasonal water use budgets varied greatly from soil to soil. These differences were
attributed to differences in plant size. Values for total monthly crop water use and seasonal
use are tabulated in Table 4.2. Cumulative use for the final eight weeks of the growing
season are shown in Figures 4.94.11
Table 4.2. Monthly and Seasonal Water Use.
May June Seasonal
Soil Type
(mm) (L/m) (mm) (L/m) (mm) (L/m)
Builder's Sand 41 53 50 64 91 117
Astatula f.s. 77 99 75 96 152 195
Pomona sand 145 186 195 250 340 436
Depths (mm) CUMULATIVE WATER USE
350
30 I from tri ntlon
0 from watr table
50
02 08 15 22 29 05 12 19 26
MAY JUNE
Builder's sand
Figure 4.9. Cumulative water use for the plants growing in Builder's sand.
Depths (mm) CUMULATIVE WATER USE
350
150
02 08 15 22 29 05 12 19 26
MAY JUNE
Astatula f.s.
Figure 4.10. Cumulative water use for the plants growing in Astatula f.s.
Depths (mm) CUMULATIVE WATER USE
350
300
250
200
150
100
50
0
02 08 15 22 29 05 12 19 26
MAY JUNE
Pomona sand
Figure 4.11. Cumulative water use for the plants growing in Pomona sand.
Redistribution and Extraction Patterns
The data acquisition system proved extremely effective in determining redistribution
and extraction patterns in the soil profile. These patterns, however, were limited to the top
300 mm of the soil profile because that is where 9 of the 10 tensiometers were located
(Figures 3.5 and 3.6)
The data acquisition system read the instruments at one minute intervals. These
readings were passed on to the microcomputer where ten minute averages were computed and
saved to disk. These data, graphed for weekly periods, provided excellent means of
monitoring the changes in soil water status in the soil profile. Figure 4.12 is a graph of
readings from the three tensiometers in the Astatula f.s. at the 100 mm depth and at radial
distances 50, 150, and 250 mm from the irrigation lateral. This figure illustrates extraction,
irrigation, and redistribution during a week in late April when the plants were small.
Since the 50 mm tensiometer was closest to the irrigation lateral and the plants, it was
the first to respond to irrigation and extraction by plant roots. As a result, the extremes
exhibited by this curve are much larger than those of the other two curves. The first high
peak of the curve occurred in the late afternoon of April 27. Immediately following is a
period of overnight redistribution of soil water succeeded by the reinitiation of extraction
the next morning. These events are marked by the U shape in the curve. At the top of the
second peak, an irrigation occurred, creating the sudden increase in soil water potential.
Figure 4.13 isolates the 24hour period around the irrigation event and permits a closer view
of how the tensiometers reacted to the wetting front.
From Figure 4.13, there was a considerable time lag between the irrigation event and
the tensiometers response to the wetting front. As expected, the lag was a function of
distance. The tensiometer 50 mm from the lateral responded to the wetting front about 10
minutes after the irrigation event which occurred at 10:30 and the tensiometer 150 mm away
responded to it about 40 minutes later. It was difficult to distinguish when the wetting front
reached the 250 mm tensiometer because the tensiometer's response was only on the order of
2 kPa. However, the tensiometer did eventually react to the approaching wetting front.
APRIL
Daily Scale
Figure 4.12. Tensiometer readings versus time on a daily scale.
ASTATULA F.S.
15 0 adial dist 5 = 550 
f radio dist. =150mm
rdia dist. = 250mm
10a
Hourly Scale
Figure 4.13. Tensiometer readings versus time on an hourly scale.
Another series of extractionirrigationredistribution cycles which occurred in the
Pomona sand during the week of June 6 is shown in Figure 4.14. In this case, the plants were
mature and once extraction began, low soil water potentials were reached very quickly. The
events of June 11 and 12 were chosen for closer study because irrigation occurred late in
the afternoon, shortly before sunset. This condition allowed redistribution and extraction to
be studied independently. Redistribution proceeded overnight with minimal influence from
extraction, while extraction was studied the following morning after most of the
redistribution had been completed. All the sudden drops (increases in soil water potential)
indicate irrigation events.
The extractionredistribution cycle of June 11 and 12 was studied by using transducer
readings over a 24hour period to graph soil water equipotential lines in the bed profile.
Equipotential lines were limited to the bed because, as can be seen in Figure 4.15, 9 of 10
tensiometers were located there. This, however, was not a serious limitation to this study
because most of the redistribution and extraction events occurred in the bed. Only half of
the lysimeter profile was instrumented because the irrigation lateral, which was located along
the center line of the bed, created a plane of symmetry. For the same reason, the tensiometer
at the 500 mm depth which was actually located on the opposite half of the lysimeter
(Figure 3.5) from the other nine tensiometers, was used as shown in Figure 4.15.
Figures 4.164.21 depict the change in equipotential lines over a 24hour period in
the instrumented Pomona sand lysimeter for the June 11 and 12 irrigation, redistribution, and
extraction cycle. Figure 4.16 shows the soil water status in the bed when the irrigation event
occurred at 17:00. Low potentials were found mostly in the top left quadrant of the profile
where, as expected, most of the roots were concentrated. The tensiometers on the right hand
side of the bed were in relatively moist soil at 10 kPa.
The next figure (4.17) shows the equipotential lines 20 minutes after the irrigation
event. By that time, the wetting front had moved past the shallowest tensiometer 50 mm
from the lateral and was beginning to affect the shallowest tensiometer at the 150 mm radial
distance. The front had not yet penetrated to the 200 mm depth.
POMONA SAND
60
S radial dist. = 50mm
S 1 t radial dist. = 150mm
S50 A radial dist. = 250mm
40
c
60
0 30
0
C 20 I A
* 10
0  __)
Daily Scale
Figure 4.14. Tensiometer readings versus time on a daily scale for the
week of June 6.
O 50 150 250 Radial Distances
and Depths (mm)
00
00 m m m TENSIOMETER
LOCATIONS
0O .a
WATER
TABLE
300 500 750
Figure 4.15. Side view of a lysimeter showing the bed, the water
table, and the tensiometer locations.
Figure 4.18 shows the equipotential lines at 18:00, one hour after irrigation. By that
time, the wetting front had clearly reached the 6 tensiometers closest to the plane of
symmetry. The equipotential lines in this figure appear to be in some disarray. This was
caused by assuming that this line source problem was twodimensional. In fact, it was not
completely so because the wetting front did not move uniformly through the soil profile
along the whole length of the bed. Since the tensiometers were not located in a single plane
as shown in the figures, but were located along the length of the bed, some instruments
responded to the wetting front earlier or later than they would have in a true 2dimensional
homogeneous, isotropic soil. This situation occurred only shortly after an irrigation event
when most of the instruments were affected by the moving wetting front.
Figure 4.19 shows the equipotential lines following overnight redistribution. This
figure distinctly shows which sections of the soil profile were affected by the irrigation
event. The potentials in the top right quadrant of the bed remained unchanged despite
irrigation while the remainder of the bed returned to field capacity. Clearly, the radial
movement of the wetting front was limited to less than 250 mm. The vertical movement of
the front, on the other hand, was more than 300 mm.
Figures 4.20 and 4.21 show the effects of extraction. At 12:00, the effects of
extraction are apparent. At 17:00, 24 hours after the irrigation event, soil water status in the
profile was very similar to the conditions that existed before irrigation. What is also apparent
from these figures is that there was very little extraction taking place near the edges of the
bed.
Overall, very little change in soil water status took place at radial distances of more
than 200 mm from the irrigation lateral. On the other hand, redistribution and extraction
were evident for at least 300 mm in the vertical direction. Similar redistribution and
extraction patterns were obtained for the builder's sand and the Astatula f.s. Radial
movement of the wetting fronts in these two soils, however, was less than in the Pomona
sand. These observations were a function of the hydraulic properties of the sandy soils
studied.
0 50 100 150 200 250 300 350 400 450 5
Distance from lateral (mm)
Figure 4.16. Equipotential lines in a Pomona sand profile at irrigation (kPa).
redistribution
17:20
0 50 100 150 200 250 300 350 400 450 500
Distance from lateral (mm)
Figure 4.17. Equipotential lines in a Pomona sand profile during redistribution (kPa).
redistribution
18:00
50
E 100
150
S200/
250
50 _5 Pomona sand
300 I
0 50 100 15L 200 25 300 350 400 450 500
Distance from lateral (mm)
Figure 4.18. Equipotential lines in a Pomona sand profile during redistribution (kPa).
3 redistribution
07:00
50
100 \
E * \
E
150 \
CI\
S200
250.
250 / Pomono sand
300 *  
0 50 100 150 200 250 300 350 400 450 500
Distance from lateral (mm)
Figure 4.19. Equipotential lines in a Pomona sand profile following redistribution (kPa).
extraction
12:00
4\
Pomona sand
0 50 100 150 200 250 300 350 400 450 500
Distance from lateral (mm)
Figure 4.20. Equipotential lines in a Pomona sand profile during extraction (kPa).
50
100
E
150
c
c 200
 I I I I I I I I I I I
0 50 100 150 200 250 500 350 400 450 500
Distance from lateral (mm)
Figure 4.21. Equipotential lines in a Pomona sand profile during extraction (kPa).
Evaporation from the Soil
Evaporation from the bare soil (furrows) between the mulched beds of each lysimeter
could not be calculated from ET data collected during the growing season. A short
experiment was therefore conducted after the growing season to determine the magnitude of
the evaporation term.
The plants were removed from the lysimeters and perforations in the plastic mulch
were sealed to eliminate transpiration or evaporation losses from the bed. The water table
was continuously recharged by the float valves, and the water table tanks were replenished
after two or three day measurement periods.
The depths of water lost to evaporation from each lysimeter was calculated by
dividing the volumes required to refill the water table tanks by the furrow surface areas of
the lysimeters (0.255 m2). The evaporation depths were then compared to Penman ET values.
Figure 4.22 presents the data for approximately three weeks in October and November, 1988.
The data are given as depths of water lost to evapotranspiration as computed by the Penman
equation and depths of water lost to evaporation from the three soils used in the study during
the two or three day measurement periods. The evaporation depths of all three soils were
generally an order of magnitude smaller than the corresponding ET values but followed the
same trend. The evaporation values were the largest for Pomona sand, followed by
Astatula f.s., and builder's sand. Consequently, it was concluded that the evaporation values
were a function of the soil properties, and in particular the particle size distribution and
organic matter content. Because the two replications of each soil produced almost identical
data, only the average data are presented in Figure 4.22.
Soil evaporation coefficients, C,, were developed for each of the soils. The
coefficient of each soil was computed by dividing the observed evaporation by Penman ET
for each of the nine measurement periods discussed above and averaging the nine values to
obtain a mean. The coefficients were C, = 0.037, 0.070, and 0.097, for builder's sand,
Astatula f.s., and Pomona sand, respectively. Daily bare soil evaporation from the furrows
Penman ET
20
E
. 10
5
12 14 17 19 21 24 26 28 31 04
Oct. Nov
Astatula f.s.
20
E 15
E
o to
S5.
12 14 17 19 21 24 26 25 31 04
Oct. Nov.
Measurement Periods
Builder's sand
20
15
10
5
0 F,
12 14 17 19 21 24 26 28 31 04
Oct. Nov.
Pomona sand
20
12 14 17 19 21 24 26 28 31 04
Oct. Nov
Measurement Periods
Figure 4.22. Comparison of Penman ET and observed bare soil evaporation
from the lysimeters.
68
during the growing season was approximated by multiplying the corresponding daily Penman
ET values by the appropriate coefficients. Figure 4.23 presents daily Penman ET pan values
and estimated bare soil evaporation values from the lysimeters for each of the three soils
during the last two months of the growing season.
800
1 15 22 29 5 1:I 19 26
MAI'
Figure 4.23. Comparison of daily Penman ET and estimated daily bare
soil evaporation from the lysimeters during the final two
months of the growing season.
'T E
A "
CHAPTER 5
NUMERICAL SOLUTION
The Richards equation was solved with a finite difference approximation. Because
line sources produce a continuous wetted band along the length of the lateral (the yaxis),
they can be studied as twodimensional problems in the xz plane. A twodimensional
numerical model was therefore developed to simulate infiltration, redistribution, and
extraction of water in the soil profile.
Several assumptions were employed in the development of the numerical model to
simplify the complexities associated with mathematically describing water movement in
unsaturated porous media. These assumptions were:
1. Physical properties of the soil were homogeneous, isotropic, and nondeforming
under varying moisture content and pressure head.
2. Physical and chemical properties of the soil were constant in time.
3. Hydraulic conductivity, K(O), and pressure head, 0, were singlevalued
functions of water content, 0, (nonhysteretic).
OneDimensional Model Development
A onedimensional vertical simulation model was developed as the precursor of the
twodimensional model. Its purpose was to verify that the finite difference equations closely
approximated the governing partial differential equation. For onedimensional flow, equation
(2.19) reduces to
at a r f 1(@)
C()= K() 1 S() (5.1)
at az ie iaz
where z = the vertical dimension.
70
Equation (5.1) was discretized with an implicit finite difference scheme that used
explicit linearization of the soil parameters as suggested by Haverkamp et al. (1977).
Application of equation (5.1) to the soil was accomplished by superimposing a mesh over the
soil profile with a finite number of grids or cells as shown in Figure 5.1. For simplicity, the
parameters K(V), C(O), and S(4) were represented by K, C, and S, respectively. The
superscript j was used as the space increment and Az as the grid spacing in space.
Depths q.urf
j=1 soil surface
j=l Az. 01 K, C,
jI *. t^j.1 Oj.1 Kj_j Cj_
K C
j=j+l Kj+ 9, + Cjp
j=j *j j Kj C
j=J 3 B K1 Ct
j=J+ water table
j=J+l bobtton=m 0=o K=K, C=0
2
Figure 5.1. The finite difference grid system used for the onedimensional model.
Adapted from Stone (1987), p.38.
The pressure head and soil parameters of each cell were computed and assigned to
the node located at the center of each cell. The total number of interior nodes was J with
the water table, which acted as the lower boundary, located at j=J++. The exterior node, at
j=J+1, was used in the finite difference equation that described the lower boundary
condition.
The numerical form of equation (5.1) was obtained by letting
a [ (a) Ow
 [K() I 1I (5.2)
Oz L Oz az
and discretizing the right hand side of equation (5.2) around j to obtain
I J Wj+ W (5.3)
'9Z i Az
Equation (5.3) was then substituted into equation (5.2) to give
1 fr 1 1 F (a4 1 1
K]j [ 1 K I I (5.4)
AzL J J L z J
The derivatives of i with respect to z were approximated as
' = 1 (5.5)
az j+i Az
ao 1 'JiI
1 1 (5.6)
Ir J Az
which when substituted back into equation (5.4) gave
A'z [K [  z [ 1 (5.7)
Az A~z Az
The superscript n was used as the time coordinate and At as the grid spacing in time.
The implicit formulation of equation (5.7) was obtained by evaluating at the n+l time level.
The time derivative in equation (5.1) was then replaced by a backwarddifference
approximation relative to the n+l time level. By using explicit linearization, the soil
parameters are known at each time step and are evaluated using the values of 0 at the n time
step (tn).
on f[ +[ 1
C p K *
At Az I Az J
[K  1 } (5.8)
The hydraulic conductivity K(tti) was taken as the arithmetic mean for all interior
nodes I s j s J where
K  + (5.9)
2
KI + K.
K = 1 (5.10)
2
Since explicit linearization was used to evaluate the soil parameters, the only
unknowns in equation (5.8) were terms containing "n+'. The equation was then rearranged
so that the unknowns were on the left hand side and the knowns on the right hand side
l j +i+ [C +(K 1 + K? + n
Az J LAt Az' J AZ +
C';I (KjI K +)
= Sn (5.11)
At Az
or in a more usable form
Aj+' + Dj + BjO: = R (5.12)
where
K"
A = (5.13)
Az'
At
K+
B = I (5.15)
AZ2
R (K_ K ) (5.16)
At Az
Equation (5.12) was written for each of the interior nodes 2 j < JI in the soil
profile. The top and bottom grids of the system required equation (5.11) to be modified so
that it would describe the boundary conditions.
Boundary Conditions
Two types of boundary conditions are usually associated with the flow equation in
a soil profile. A Dirichlet condition specifies the dependent variable (in this case the pressure
head, 4). A Neuman condition specifies the derivative of the dependent variable which for
the soil water problem means a specification of the flow through the boundary (Feddes et al.,
1978). To simulate the conditions encountered in this study both types of boundary
conditions were specified.
Too boundary. Physically, the top boundary of the soil profile consisted of the plastic
mulch. The only flow through the soil surface allowed by the mulch was surface infiltration
due to irrigation. Mathematically, this was modeled by a flux boundary condition with a
surface flux, qurfqe (q,). The surface flux during irrigation was set equal to the water
application rate until the irrigation event was completed. At all other times, q, was set to
zero. The total depth of water infiltrated at any time was calculated as the summation of the
incremental infiltration volumes at all previous time steps. A flux across a boundary in the
vertical direction is described by the partial differential equation
q=K (5.17)
L8z J
When equation (5.16) was discretized around j it became
qj Kj l1 (5.18)
k Az
which when evaluated at j = 1, the soil surface (Figure 5.1), gave
I = K I q (5.19)
L Az
When equation (5.7) was evaluated at the first node, j=l, it became
n1 rrn f I On+1
C  K~,2 I
At Az Az
[K/2 I  k Sn (5.20)
L Az
The third term of equation (5.20) is the same as equation (5.19). Substitution of q, for the
term reduced equation (5.20) to
Cn = K/2 1 +q, S' (5.21)
At Az [ I Az
When the equation (5.21) was rearranged so that unknowns were on the left and knowns on
the right it became
[Cn K/2"I K n1 K1 n
I i S; (5.22)
At Az2 A At Az
or in a more usable form
Dl 'b+ + Bln+1 = R (5.23)
where
c;
D, B, (5.24)
At
Kn/
B, = (5.25)
Az2
CnV q, Kn/.
RI=  S (5.26)
At Az
Bottom boundary. A constant potential boundary condition was used to simulate the
water table. When equation (5.11) was evaluated for the last node, j = J, it became
[ K_ 1 + (Kj + Kb+ p) + K In 1l
AZ2 t AZ2 Z2 J
Cn n (Kjn_ Knj,)
S; (5.27)
At Az
Because J+1 was the exterior node at the lower boundary which was held at a constant
pressure head, ,botom,, "+ = .+1 = bottom. Therefore, the term containing Int in
equation (5.27) was known and could be moved to the right hand side of the equation, which
when simplified, became
76
AjOb. + Dj;+1 = Rj (5.28)
where
A = (5.29)
Az2
Cn Kn
Dj= Aj + (5.30)
At Az'
Cnn, (Kn Kn,) K+,,
R = + bottom S; (5.31)
At Az AzI
The flux across the lower boundary, bottom (qb) was calculated by evaluating
equation (5.18) at the lower boundary (j=J+i). Physically, qb represented the rate at which
water moved from the water table into the soil profile or drained to the sump.
An equation for a flux boundary condition at the lower boundary was also developed
to allow the model to be compared to solutions developed by other researchers. The equation
was derived by following the same steps used to develop equation (5.22).
Solution
The system of equations which were formed by applying equations (5.28), (5.12), and
(5.23) to the appropriate nodes produced a tridiagonal matrix equation which was implicit in
terms of pressure head, '.
D, B 1 R,
A2 D, B, 22 R,
= (5.32)
A_J Dj_ Bj jj Rj
A, Dj 0 Rj
77
The tridiagonal matrix was readily solved by a Gaussian elimination method. In this
work, a very concise algorithm utilizing recursion formulas for the solution of the tridiagonal
system which was presented by Carnahan et al. (1969), was used.
A computer program used the finite difference equations and solved the tridiagonal
matrix to simulate infiltration, redistribution, extraction, and upward flow from a water
table. The program was written in Microsoft FORTRAN and designed to run on
microcomputers.
Time Step
Implicit numerical techniques generally permit the use of much larger time steps
than explicit techniques while maintaining stability. A suitable time step may be determined
by trial and error and kept constant throughout the simulation. A much more efficient
approach, however, is to have a variable time step which decreases in size during periods of
increased soil water movement (infiltration) and increases during periods of relative inactivity
(redistribution). Zaradny (1978) proposed that the time step be estimated by
( Az
At < (5.33)
1qmaxl
where
At = the variable time step (hrs),
Az = grid spacing in space (mm),
lmqx = the maximum net flux into any grid (mm/hr), and
r = the maximum permissible change in water content, 0, for any grid in the
soil profile where 0.015 < f < 0.035.
Feddes et al. (1978) used the lower value of f to accommodate the rapid movement of water
during infiltration and early stages of redistribution. Higher values of { were used during
slow changes in boundary conditions such as continuous upward flow of water or during the
78
later stages of redistribution. Clark (1982) and Stone (1987) also used this method to compute
time steps.
Updating of Soil Parameters
The soil parameters 0, C, and K corresponding to the new values of 0 at the end of
every time step were determined from equations (2.6), (2.8), and (2.9), respectively, for the
three soils used in the lysimeters (builder's sand, Astatula f.s., and Pomona sand). Because
the observed soilmoisture retention curves and those predicted by equation (2.6) were well
matched (Figure 3.11), this approach was much more efficient and about as accurate as a
table lookup procedure.
The tabulated parameters of five other soils presented by Clark (1982) were included
in the model to give it more flexibility. A tabulated data search was used to determine the
corresponding values of K, C, and 6 for these soils. Because the search had to be conducted
for each grid at the end of every time step, the logarithmic interpolation method used by
Clark (1982) was employed to reduce computer time. The table was set up with linearly
increasing multiples of the logarithmic value of 0. Since the tabulated soil parameters were
discrete values, errors arose from the interpolation procedure. The errors however, proved
to be less than 1% and were neglected. The tabulated parameters of the five soils, (Arredondo
f.s., Rehevot sand, Adelanto loam, Pachappa loam, and Yolo light clay), are presented in
Appendix C.
Mass Balance
A mass balance was computed at the end of every time step and used to check the
stability of the simulation model. The water stored in each grid was integrated throughout
the profile and compared to the cumulative inflows and outflows. Inflows consisted of
surface infiltrations and upward flow from the water table. Outflows included drainage,
root extractions, and evaporation from the soil surface which was a component of the surface
grids' root extraction term. The mass balance error was calculated as the difference between
the water stored in the profile at any time and the sum of inflows and outflows.
OneDimensional Model Verification
The accuracy of the onedimensional model was determined by comparing it to
experimental data and analytical and numerical solutions found in the literature. The
comparisons were selected to test every aspect of the model with soils that had widely varying
hydraulic properties.
Analytical solutions of the Richards equation generally deal with steadystate
situations or simplified boundary conditions or geometries, but are useful because they
provide a true solution against which numerical simulations can be compared. Gardner (1958)
developed several analytical solutions based on the matric flux potential concept (Philip, 1959)
to describe steadystate evaporation from a water table in soils whose hydraulic conductivity
function could be described by
a
K(O) = (5.34)
(, + b)
where a, n, and b are constants.
Gardner and Fireman (1958) used one of these analytical solutions to calculate steady
state evaporation from a Yolo light clay with a water table 1050 mm from the soil surface.
To obtain the solution, they empirically fitted equation (5.34) to the Yolo light clay hydraulic
conductivity data that Moore (1939) presented in his landmark paper. From the resulting
curve they determined that a = 400 cm3/day, b = 400 cm2, and n = 2 and used these constants
in the analytical solution that they presented. They simulated steadystate evaporation at the
rate of 0.033 mm/hr (0.08 cm/day) and compared their results to the data given by Moore
(1939) for steadystate evaporation from a column of Yolo light clay with a water table at a
depth of 1050 mm.
The onedimensional model developed in this study was used to simulate steadystate
evaporation for the same conditions. For this validation, the evaporation rate was simulated
as a negative surface flux and K(i) was computed from equation (5.34). The other soil
80
parameters were updated with the tabulated data search. The numerical solution was in
excellent agreement with the analytical solution and the experimental data as shown in
Figure 5.2.
Simulated infiltration into an initially very dry soil profile was compared to a study
by Rubin and Steinhardt (1963). They developed a model to study constant intensity rainfall
infiltration into Rehevot sand by implicitly solving the 0form of the Richards equation.
They presented results for rainfall intensities of 12.7 and 47 mm/hr on a semiinfinite soil
profile with a uniform initial soilmoisture content of 0.005.
The onedimensional model developed in this work was used to simulate infiltration
and redistribution for identical initial and boundary conditions. Soil parameters were updated
with the tabulated data search using the Rehevot sand soil hydraulic characteristics table
(Appendix C). Results were in excellent agreement with Rubin and Steinhardt (1963) and are
presented in Figures 5.3 and 5.4 in terms of soil water content profiles.
50
300
800 su  This work
 Gardner and Frean (1958)
SMoore (1939)
1050
0 2000 4000 6000 8000
Pressure Head, i, (mm)
Figure 5.2. Simulated, analytical, and experimental results of steadystate evaporation
at a rate of 0.033 mm/hr from a water table in Yolo light clay.
50
1.56 hrs
150
250 hrs
250
350 4.38 hrs
450
6.25 hrs
550
650  This work
r Rubin and Steinhardt (1963)
750 
0.00 0.05 0.10 0.15 0.20 0.25 0.3'
Soil Moisture Content, 0
Figure 5.3. Simulated results of soil water content profiles for infiltration
into Rehevot sand under constant rain intensity of 12.7 mm/hr.
S50
150 0.66 hrs
150
1.00 hrs
S 250
E 350 66 hrs _
S 450 2.33 hrs ,
550
550  This work
750
0.00 0.05 0.10 0.15 0.20 0.25 0.
Soil Moisture Content, 6
Figure 5.4. Simulated results of soil water content profiles for infiltration
into Rehevot sand under constant rain intensity of 47 mm/hr.
82
Warrick (1974) also solved the onedimensional, steadystate moisture flow equation
using the matric flux potential to obtain a linearized form. He presented solutions for a
semiinfinite flow medium and a finitedepth medium overlying a shallow water table with
a flux surface boundary condition and arbitrary plant water withdrawal functions. He
assumed unsaturated hydraulic conductivity, K(0), to be of the form
K(O) = Ke(ot) (5.35)
where K, is the saturated hydraulic conductivity and a is an empirical constant depending
on soil texture.
Warrick (1974) used an arbitrary soil with K, 4.17 mm/hr and a = 0.001 mm' to
present numerical examples of his analytical solutions. One set of numerical results was
developed for a soil profile with a water table 1000 mm from the soil surface. The arbitrary
extraction function was set to 0 and surface fluxes of 0.0 mm/hr and 1.5 mm/hr were used.
Equation (5.35) and a linear soilmoisture release function were incorporated into the
onedimensional model to allow comparison with Warrick's solutions. Steadystate pressure
head were simulated for identical boundary conditions and they were in excellent agreement
with the analytical solutions. Both solutions and are shown in Figure 5.5.
In order to evaluate the performance of extraction in the model, the simplest water
extraction function given by Warrick (1974) was chosen and coded into the model. This was
a constant sink term give by
S(z) = a (5.36)
where a is the constant extraction rate per unit length. The total soil profile water extraction
rate used by Warrick (1974) was 0.42 mm/hr (S=aL where L is the depth to the water table).
The numerical simulations were run to steadystate using the same conductivity
function and initial and bottom boundary conditions as before but with infiltration rates of
0.42 and 0.63 mm/hr. The results were in excellent agreement with the analytical solution
and are shown in Figure 5.6.
E
E 400
q,=0.0 mm/hr
S600
This work
 Warrick (1974)
1000
0 200 400 600 800 1000
Pressure Head, (mm)
Figure 5.5. Simulated and analytical results of steadystate infiltration above
a shallow water table at rates of 0.0 and 1.5 mm/hr.
200 400 600 800 1000
Pressure Head, "1 (mm)
Figure 5.6. Simulated and analytical results of steadystate infiltration
above a shallow water table with a constant extraction rate
of 0.42 mm/hr.
q,=0.63 mm/hr
q,=0.42 mm/I
This work
Warrick (197
hr
4)
TwoDimensional Model Development
The process of infiltration, redistribution, and extraction under line sources or closely
spaced emitters can be treated as a twodimensional problem. The governing equation for
twodimensional flow in a porous medium is
C()  K() + [K(0) I S(') (5.37)
at ax L ax az I
where x represents the horizontal dimension along the xaxis and z represents the vertical
dimension along the zaxis.
As shown in Figure 5.7, the line source lies on the soil surface along the yaxis and
intercepts the xz plane at the origin. The trickle discharge per unit length of the line source
is given by Q(t) in mms/(mmhr).
Equation (5.37) was solved numerically with a Gaussian elimination implicit finite
difference scheme and explicit linearization of soil parameters. The flow domain was divided
into a mesh of rectangular cells having dimensions Ax, Az, and unity in the y direction. The
index i (i=l, 2, ..., I) was used to locate the cells in the xdirection with respect to the line
source, and j (j=l, 2,.... J) was used to number the cells in the vertical direction (Figure 5.7).
Horizontal flow in equation (5.37) is described by
K() (5.38)
ax L ax J J ax
The finite difference approximation of this equation was obtained by discretizing the right
hand side of equation (5.38) around i to obtain
i Wi+.j W i (5.39)
Lax ij Ax
Equation (5.39) was then substituted into equation (5.38) to give
yl sou
line source
line of
symmetry
soil surface
Figure 5.7. Schematic diagram of the finitedifference grid system for the
twodimensional model of water movement and extraction.
*s,>
ij1
il'j ij i+lj
i,j+l
K. [ [K i (5.40)
The derivatives of 0 with respect to x were approximated as
S J iJ (5.41)
ax i+j Ax
] ij 'ij (5.42)
ax ij Ax
which when substituted back into equation (5.40) gave
^(r j11 r r '"111
1[Ki+u [ A ii K iJ x J 1 (543)
Equation (5.43) represents the finite difference approximation of horizontal flow
through a porous medium. This equation was combined with the finite difference
approximation of vertical flow previously derived for the one dimensional model (equation
(5.7)) to produce the twodimensional finite difference equation for flow through a porous
medium.
S on
CAx
At z [ Axz
[K I' " 1
j S1j (5.44)
LKJ A J] 5z
The equation was then rearranged so that the unknowns were on the left hand side and the
knowns were on the right hand side.
[K ~t Kn ,
Ax Az/
C ij Ktfj + K KP+ K
+ + +  [ + i P+
L t Ax2 Az2~
oi+,lj 
[Ax2 AZ
c? K Kj+
= J + Si7
At Az
or in a more usable form
A + Bi,j ',j + Diji + E i,j + Eijj+J i = ij
where
K?,j
Ai,j = 
AZ2
K"
Ax2
C"
Di,j = Ai, Bi E,j Fj
At
SK"
Ei,j = 
Ax2
(5.45)
(5.46)
(5.47)
(5.48)
(5.49)
