MODELING THE EVAPORATION AND TEMPERATURE
DISTRIBUTION OF A SOIL PROFILE
BY
CHRISTOPHER LLOYD BUTTS
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
1988
112 S
ACKNOWLEDGEMENTS
I would like to acknowledge, first and foremost, GOD's blessings
in giving me the talents and insights with which to complete this
dissertation. My wife, Sherry, has provided steadfast support and
encouragement without which this task would have been impossible. I
also thank my committee chairman, Dr. Wayne Mishoe, for his guidance
and encouragement and the free reign to conduct the research as I
deemed necessary. I thank my committee, Drs. James Jones, Khe Chau,
Hartwell Allen, Calvin Oliver and Mr. Jerome Gaffney, for their
guidance and support. To the Agricultural Engineering technical staff
in the fabrication and instrumentation of the weighing lysimeters
constructed for this project, I offer my thanks. I thank Bob Bush for
his diligent maintenance and operation of the weighing lysimeters.
Finally, I acknowledge the contribution to my wellrounded education of
my fellow graduate students, particularly, Bob Romero, Ken Stone, Matt
Smith, Ashim Datta and Kumar Nagarajan.
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS . . . . . . . . . . . . ii
KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . v
ABSTRACT . . . ...... . ....... . . . . . .x
CHAPTERS
I. INTRODUCTION . . . . . . . . . . . 1
Problem Statement . . . . . . . . 1
Research goals and objectives . . . . . . 3
II. ENERGY AND WATER TRANSFER IN A SANDY SOIL: MODEL DEVELOPMENT 5
Introduction . . . . . .
Literature Review . . . .
Model Objectives . . . . .
Model development . . . .
Determination of Model Parameters
Numerical solution . . . .
III. EQUIPMENT AND PROCEDURES FOR MEASURING THE MASS AND ENERGY
TRANSFER PROCESSES IN THE SOIL .
Introduction . . . . .
Objectives . . . . . .
Lysimeter Design, Installation a
and Construction . . .
Experimental Procedure . . .
Data Analysis . . . .
Experimental Results . . .
Summary . . . . . .
. . . . . . .o
. . . . . . .
d Calibration Design
. . . . . . .
. . . . ... .
. . . . . . .
. . . .
. . . . .
IV. MEASURING THERMAL DIFFUSIVITY OF SOILS . .. . . .
Introduction . . . . . . .. . . .
Objectives . . . . . . .. ....
Literature Review . . . . . . . . .
Procedure . . . . . . . . . . .
Results and Discussion . . . . . . ...
Conclusions . . . . . . . . . . .
5
6
19
20
30
43
66
66
67
67
85
89
90
96
134
134
135
135
141
143
148
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
ooooooooo
oeooeooooo
oooooooooo
oooeoooooo
eoooooooo
oooooooeeo
V. MODEL ANALYSIS . . . . . . . . .. . .
158
Introduction . . . . . . . . .. . 158
Validation . . . . . ...158
Sensitivity Analysis . ..... . . . .. 170
Summary . . . . . . . . . . 176
VI. SUMMARY AND CONCLUSIONS . . . . . . . . . 208
BIBLIOGRAPHY . . . . . . . . .. . . . . 212
APPENDIX A ADI FORMULATION OF COUPLED HEAT AND MASS
TRANSFER MODEL . . . . . . . . . . . 219
APPENDIX B METEOROLOGICAL DATA FOR LYSIMETER EVAPORATION
STUDIES . . . . . . .... . . . . . 225
BIOGRAPHICAL SKETCH . . . . . . . . . . . 241
KEY TO SYMBOLS AND ABBREVIATIONS
Variable Definition
B empirical function for determination of Cer and Chr
dimensionlesss]
Cdr surface drag coefficient dimensionlesss]
Cer Dalton number, dimensionless mass transfer
coefficient
Chr Stanton number, dimensionless heat transfer
coefficient
Cpa specific heat of moist air [J'kg1"K1]
Cps specific heat of solid constituent of soil mixture
[Jkg1.oK]
Cpw specific heat of water [Jkg1'K1]
Cpv specific heat of water vapor [Jkg1K1]
Cs volumetric heat capacity of soil mixture (solid,
liquid and air phases) [J'kgl1(1]
Da diffusivity of water vapor in air [m2"s1]
DL unsaturated hydraulic diffusivity of soil [m2.s1]
Dt thermal diffusivity [m2.s1]
Dv diffusivity of water vapor in soil [m2.s1]
dwj height of volume of soil associated with node j. [m]
dxj width of soil cell j in the x direction (unity) [m]
dyj width of soil cell j in the y direction (unity) [m]
dzj vertical distance between nodes j and j+1. [m]
E(z,t) rate of phase change of liquid water to water vapor
@in the soil as a function of space (z) and time (t)
[kg H20)m(soil)s'
Variable Definition
Ej,n rate of phase change of liquid water to vapor in
soil cell j at time step n [kg(H20)'m(oi)*s" ]
ETa actual evapotranspiration [mm]
ETp potential evapotranspiration [mm]
ea water vapor pressure in ambient air [Pa]
eas saturated water vapor pressure in ambient air [Pa]
es(T) saturated vapor pressure at temperature, T [Pa]
evs(O) vapor pressure at water potential, 0 [Pa]
G sensible heat flux rate from the soil to the air
[Wm ]
g gravitational acceleration [m's2]
gj shape factor in the jth principal axis for thermal
conductivity calculations dimensionlesss]
H vertical component of heat flux in the atmosphere
[Jm ]
hfg latent heat of vaporization of water [J3kg1]
hh convection heat transfer coefficient for soil
surface [W m'2Ki]
hm convection vapor transfer coefficient for soil
surface [ms ]
K hydraulic conductivity [m0'm2.s1]
Kh, Km, Kw turbulent transfer coefficients for peat, momentum
and water vapor, respectively [m s ]
kij ratio of temperature gradient in ith soil
constituent to the temperature gradient in the
continuous constituent (water or air) in the
direction of the jth principal axis dimensionlesss]
LE latent energy transfer in the atmosphere [Jm2]
n number of moles of gas present, used in ideal gas
law [mol]
Variable Definition
P precipitatior: flux of water at soil surface
[m H20)m s]
p atmospheric pressure [Pa]
Po reference atmospheric pressure [Pa]
Qv vertical flux density of water vapor in the
atmosphere [kgm 's ]
q specific humidity of ambient air [kgapor/kgdry air]
qcj+l,n rate of heat flux conducted into node j from node
j+1 at time step n [Wm ]
qej,n rate of latent heat loss due to evaporation at node
j and time step n [W'm" ]
qLj+1,n heat flux rate into node j from node j+1 via
movement of liquid at time step n [W'm2]
qsj,n rate of sensible heat change of node j at time step
n [W'm3]
qvj+1,n rate of heat carried into node j from node j+1 by
water vapor movement at time step n [W'm2]
R universal gas constant [kgm2s2.mol1(o1]
Rw gas constant for water vapor, determined by dividing
the universal gas constant, R, by the molecular
weight of water [m2.s2OK1]
Re Reynolds number dimensionlesss]
Ri Richardson number dimensionlesss]
Rn net r diation flux incident upon soil surface
[Wm ]
S soil porosity [m3m3]
s slope of the saturated vapor pressure line
[kPa.oK1]
Ta ambient air temperature [OK]
Tj,n soil temperature at node j and time step n [OK]
Tavgj (Tj+l,n + Tj,n)/2 [OK]
Variable Definition (continued)
t time [s]
U wind velocity [msl]
U* shear velocity [ms1]
V volume [m3]
WS Wind speed [mss1]
xi volume fraction of the ith soil component [m3m3]
z depth below soil surface [m]
z0 surface roughness length [m]
a soil tortuosity dimensionlesss]
I psychrometric constant [kPaCK1]
7r diabetic influence function dimensionlesss]
soil water potential [m]
X thermal conductivity of soil composite at node j
[Wm1.OK1]
v kinematic viscosity [m2.s1]
Pa density of moist air [kg'm3]
Ps dry bulk density of soil [kg'm3]
Pv water vapor concentration [kgH20'ms3il]
Pvj,n water va or concentration at node j and time step n
[kgH20'msoil]
Pva water va or concentration of ambient air
[kgH20'mairl
Pvs water vapor concentration in the air at saturation
[kgH20mlr]
Pw density of water [kgm3]
9 dimensionless water content
viii
Variable Definition (continued)
0 volumetric water content of soil [m3.m3]
Ojn volumetric water content of soil at node j and at
time step n [mn3m3]
Or residual volumetric water content [m3.m3]
es volumetric water content at saturation [m3m3]
7 shear stress [Nm2]
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
MODELING THE EVAPORATION AND TEMPERATURE
DISTRIBUTION OF A SOIL PROFILE
By
Christopher Lloyd Butts
December, 1988
Chairman: J. Wayne Mishoe
Cochairman: James W. Jones
Major Department: Agricultural Engineering
Evaporation of water from the soil is controlled by the transport
phenomena of energy and mass transfer. Most procedures for estimating
the loss of water from the soil assume that the soil is an isothermal
medium with evaporation of water occurring at the soil surface.
Estimating evaporation as a surface phenomenon independent of soil
hydraulic and thermal properties can lead to overestimation of the
amount of water lost from the soil.
Weighing lysimeters measuring 2 x 3 x 1.3 m3 were constructed
capable of detecting a change in weight equivalent to 0.02 mm of
water. Recorded data consisted of net radiation, soil temperature,
water content, load cell output, ambient air temperature, relative
humidity, windspeed and precipitation. Measured soil properties
included dry bulk density, porosity and thermal diffusivity.
A model describing the continuous distribution of heat, water and
water vapor within the soil was developed. The evaporation of water
was allowed to occur throughout the soil profile in response to the
assumed equilibrium between the liquid and vapor phases. Surface heat
and mass transfer coefficients were determined using equations based
upon equations of motion in the atmosphere. Thermal diffusivity of the
soil was measured and incorporated into the simulation.
Simulation results included temporal values of the cumulative and
hourly water loss from the soil, and spatial distribution of
temperature, water and water vapor in the soil. Experimental data
obtained from the weighing lysimeters were used to validate the model.
Prior to model calibration simulated soil temperatures were within 2 OC
of measured values and simulated cumulative evaporation was within ten
percent of measured water loss. Sensitivity analysis indicated that
calibration could be achieved by relatively small adjustments in the
values of the surface heat and mass transfer coefficients.
Experimental and simulated evaporation rates exhibited a diurnal
pattern in which maximum evaporation rates occurred approximately
midday then decreased to near zero at sunset. Under some conditions,
water continued to evaporate from the soil overnight at a rate of
approximately 0.03 mm/h while in some cases a net gain of water was
observed overnight.
CHAPTER I
INTRODUCTION
Problem Statement
The process by which water is lost from the soil has been the
object of considerable research over the years (Brutsaert, 1982).
Water is lost from the soil by runoff, drainage, evaporation and
transpiration. Evapotranspiration (ET) has been used to describe the
combined evaporative losses from the soil and plant. Potential evapo
transpiration (ETp) has been defined as the amount of water
transferred from a wet surface to vapor in the atmosphere (Penman,
1948) and has been used as an estimate of the maximum amount of water
lost to the atmosphere from a given area. Estimates of actual
evapotranspiration (ETa) are frequently made by the use of coefficients
which vary according to ground cover and the region of the country.
The need for knowing ETa on a field basis arises for a variety of
reasons including design and selection of irrigation equipment,
optimization of irrigation management, and determination of water
requirements for agricultural and municipal entities (Heermann, 1986).
Regional estimates of ETa are necessary for hydrologic simulations
(Saxton, 1986).
Mathematical models have been developed to provide estimates of
both ETa and ETp. The use of evapotranspiration models has expanded in
recent years due to the advent of crop growth simulations requiring the
distribution of soil water as a component of the model. One such model
used to integrate the soil water and the plant water requirements was
1
2
developed by Zur and Jones (1981) and focused on estimating water
relations for crop growth. A soil water model which simulates ETa is
an integral part of the soybean production model, SOYGRO, developed by
Wilkerson et al. (1983). One of the primary uses of SOYGRO is to
determine the yield response of soybean to various environmental
conditions such as drought.
Even though ET models are currently used extensively, many
research needs remain. At best, crop coefficients are rough estimates
computed from previous experiments with the crop of concern. There is
a need to develop a better understanding of the variability of soil,
crop and environmental factors that determine crop coefficients.
Scientists need improved procedures for determining crop coefficients
for various cultural practices such as mulching, notill and row
spacing.
One of the most pressing needs in the area of evapotranspiration
research is in the area of sensors which will detect soil water status
for a particular field (Heermann, 1986; Saxton, 1986). One of the
rising technologies for determining soil water status is the use of
remote sensing of soil temperature for determining soil water status
(Soer, 1980; Shih et al., 1986).
Evaporation from the soil is driven by the energy balance of the
surrounding environment. Because of the release of latent energy is
coupled with the sensible energy balance, it is necessary to consider
both components of the energy balance to completely describe the
evaporation process. Soil temperature has been found to influence many
processes in the soil. For example, soil temperature influenced
nodulation of roots of Phaseolus vulgaris L. (Small et al., 1968) and
3
root growth for various stages of development of soybean (Brouwer,
1964; Brouwer and Hoogland, 1964; Brouwer and Kleinendorst, 1967).
Blankenship et al. (1984) have shown that excessive temperatures in the
soil zone in which peanuts are produced (geocarposphere) can reduce
yields and increase susceptibility to disease, fungus and physical
damage. Unpublished data collected by Sanders (1988) for peanuts grown
in lysimeters with controlled soil temperatures indicated that the
maturation rate of peanuts were affected by soil temperature. The fact
that movement of water vapor in the soil in response to temperature
gradients can be significant was determined by Matthes and Bowen (1968)
and more recently by Prat (1986). Neglecting soil temperature in a
mass transfer model in the soil could cause significant errors in the
estimation of ETa from the soil if soil temperature were neglected in
the modeling analysis.
Research Goals and Objectives
An understanding of the processes of heat and mass transfer in the
soil is of interest to a variety of agricultural scientists,
particularly those involved in modeling agricultural systems and
developing management strategies. Evapotranspiration is a complex
process involving energy and mass transport in the plant and the soil
as well as the chemical processes such as photosynthesis and
respiration which require water from the soil and release water from
the reaction. It has been suggested by some researchers that the two
processes involving the plant and the soil be studied independently so
that greater insight might be obtained into each. Therefore, the goal
of this research is to study the process of evaporation from the soil
independent of transpiration for the purposes of furthering the
knowledge of the parameters affecting evaporation from the soil and
guiding future research. The specific objectives are:
1. to develop a model to simulate the evaporation of water from
the soil which couples the energy and mass transfer
processes;
2. to monitor evaporation of water and vertical distribution of
temperature and water in a sandy soil;
3. to determine the validity of the model by comparison of
simulated evaporation rates, soil temperatures and
volumetric water content to those measured in a sandy soil;
4. to determine the sensitivity of the evaporation process to
changes in various soil properties and environmental
conditions; and
5. to calibrate the model for local conditions using data
collected in Objective 2 utilizing information gained from
the sensitivity analysis.
CHAPTER II
ENERGY AND WATER TRANSFER IN A SANDY SOIL: MODEL DEVELOPMENT
Introduction
Models of varying complexity may be used to simulate or estimate
the amount of water which evaporates from the soil. The complexity of
the model should be determined by the projected use of the results and
the data available as input for the model. Empirical models generally
relate one or more parameters by regression analysis to evaporation
measured under various conditions. Variables or parameters typically
used in empirical models are evaporation pan data, air temperature, and
day length or solar radiation. Monthly estimates of evaporation for
relatively large areas are typically obtained from empirical models
(Jones et al., 1984).
Resistance analog models employ the concept of the electrical
analog to mass flow, where the flux of water vapor leaving the soil
surface is expressed as a potential (vapor pressure) difference divided
by a resistance. The resistance term for the mass flow is a function
of the boundary layer of the lower atmosphere and the mass transfer
characteristics of the soil (Camillo and Gurney, 1986; Jagtap and
Jones, 1986). These resistance models have the capability of
estimating actual evaporation from the soil on a daily or hourly basis.
Evaporation can also be estimated using mechanistic models
describing the conservation of mass, momentum, and energy in the lower
atmosphere and the soil. These models are the most detailed in their
derivation and have the advantage that they can provide substantial
5
6
insight into the actual processes describing the phenomenon of
evaporation. However, considerable data may be required as input to
provide reasonable estimates of the heat and mass transfer.
Literature Review
Empirical Models
Jones et al. (1984) reviewed some of the empirical procedures to
estimate crop water use noting that all require calibration for a
particular geographical region. Pruitt (1966) reported that serious
errors can occur when using pan evaporation data to estimate crop water
use, particularly under strong advective weather conditions. The pan
should be installed properly ensuring sufficient fetch surrounding the
pan.
The Thornwaite model, as presented by Jones et al. (1984), uses
monthly averages of air temperature and day length to estimate monthly
evapotranspiration (ET). Criddle (1966) summarized the BlaneyCriddle
model which predicts the actual water use by a crop from monthly
average air temperature and percent daylight and introduces crop
consumptive use coefficients. The BlaneyCriddle model, developed for
use in arid regions, greatly overestimates evaporation for the humid
climate of Florida during the summer months prompting Shih et al.
(1977) to replace the percent daylength with monthly net radiation to
account for the increased cloud cover.
Resistance Models
Resistance models presented in the literature are of the general
form
cpaPa (e(Ts) e(Td) )
E = Hf R (21)
where:
E = evaporation rate [ms1]
Cpa = specific heat of air [Jkgl'o1]
Pa = density of dry air [kg'm3]
hfg = latent heat of vaporization [J'kg1]
= psychrometric constant [kPa*'K1]
e(T) = saturated vapor pressure [kPa] at temperature T
Ts = soil surface temperature [OK]
Td = ambient dew point temperature [OK]
R = resistance to vapor movement from the soil to air
[kgsm4]
The resistance term for the models presented by Conaway and Van Bavel
(1967), Tanner and Fuchs (1968), and Novak and Black (1985) represents
the resistance due to the laminar boundary layer. The boundary layer
resistance is a function of the wind speed, atmospheric instability,
and the surface roughness height. Their models are used to predict
evaporation from a wellwatered bare soil surface. Jagtap and Jones
(1986) and Camillo and Gurney (1986) developed resistance terms to
include the boundary layer resistance in series with the resistance of
vapor flow in the soil. The soil resistance term is included since
after the soil surface dries, the water must change to vapor in the
soil below the surface then diffuse to the soilatmosphere interface.
The soil resistance term developed by Jagtap and Jones was determined
by regression analysis as a function of cumulative evaporation, water
in the soil profile available for evaporation, and a daily running
average of the net radiation. The net radiation empirically accounted
for the heat flux into the soil, while the ratio of the cumulative
8
evaporation to the available water parameterized the thickness of dry
soil through which the water vapor diffused. The soil resistance term
used by Camillo and Gurney (1986) was a linear regression of the
difference between the saturated volumetric water content and the
actual soil water content. Both models were intended to provide daily
evaporation rates.
The resistance model estimates daily evapotranspiration relatively
well. However, the soil resistance term varies with time depending
upon the distance through which the water vapor must diffuse from below
the soil surface. One point that should be noted is that the vapor
pressure difference used in the analog models discussed is the
difference between the saturated vapor pressure at the temperature of
the soil surface and that of ambient air. In reality, evaporation
frequently occurs below the soil surface and the vapor pressure
gradient from the zone of evaporation may be substantially different
from that at the surface.
Mechanistic Models
Models derived from the basic physical relationships have the
advantage of providing estimates of evaporation on a relatively short
time scale. However, considerably more input data describing the
surface boundary conditions as well as the soil thermal and hydraulic
properties are required. Theoretical models have been developed from
four different perspectives (King, 1966; Goddard and Pruitt, 1966;
Fritschen, 1966; Penman, 1948). They are categorized as follows
1. a mass balance profile method,
2. a mass balance eddy flux method,
9
3. an energy balance method, or
4. a combination method.
The two mass balance methods (1 and 2 above) focus on the
equations of motion for the atmosphere and require accurate measurement
of velocity and temperature profiles or mass fluxes. The equations of
conservation of momentum, mass, and energy were used in the mass
balance methods to develop expressions for the rate of evaporation of
water and the shear stress between to vertical positions in the lower
atmosphere as
Qv = PaKw q (22)
az
S = p aKm aU (23)
az
H = cpaaKh aT (24)
az
where:
Qv = vertical flux density of water vapor [kgm2.s1
Pa = density of air [kg'm3]
7 = shear stress rate [N'm2]
H = vertical flux of heat [J.m2]
z = vertical distance [m]
q = specific humidity [kgapor/kgdry air]
U = velocity [m's1]
Cpa = specific heat of dry air [J*kg1Kl]
10
T = air temperature [OK]
Kh, Km, Kw = turbulent transfer coefficients for Peat, momentum
and water vapor, respectively [m' s ]
Combining the equations 22 and 23 and rearranging the following
expression for evaporation is obtained.
Q V Kwm ~ (25)
The profile methods determine the rate of evaporation by analysis of
the vertical profile of various atmospheric variables such as specific
humidity or vapor pressure and wind velocity. King (1966) described
the procedures for determining the evaporation rate for adiabatic wind
profiles (neutral stability) and stratified conditions. One of the
basic assumptions used in the profile methods is that the turbulent
transfer coefficients of water vapor and momentum are equal. This
implies that the momentum and mass displacement thicknesses are the
same. The Richardson number (Ri) or the ratio of height (z) to the
MoninObukov length (L) are used as a measure of atmospheric
instability and can be determined by equations 26 and 27,
respectively.
aT
T ( +r)
Ri = aU (26)
au
L (27)
( kgH
cpaqaT
where:
= acceleration due to gravity [ms2]
r = adiabatic lapse rate [(K(m1]
S 9.86 x 103
k = von Karman's constant dimensionlesss]
= 0.428
U* = shear velocity [mns1]
For neutral conditions the familiar logarithmic velocity profile was
used. For unstable or diabetic conditions a variation of the
logarithmic profile was used or the KEYPS function (Panofsky, 1963) was
used to determine the velocity profile. King stressed that in using
profile methods to determine evaporation, one must take extreme care in
measuring wind velocity. He suggested that spatial averaging be used
for wind speed measurements near ground level and that 30 to 60minute
averages be used considering the steady state assumptions made in
developing the equations for the profile equations. Corrections for
lower atmospheric instability caused by density gradients in the air
near the earth's surface must also be used.
Another method requiring only measurements of atmospheric
parameters is the eddy flux method (Goddard and Pruitt, 1966). The
method is based upon the turbulent transport equations of motion in
the atmosphere. The equations of turbulent motion in a fluid include
the transient random fluctuations in velocity, temperature, and mass
concentration of water vapor of the air. The following equations
describe the shear stress (T), sensible heat flux (H) and latent heat
flux (LE) due to the turbulent motion of the atmosphere
au
7 = PaKm, awu (28)
H = cpaPaKh + cpaPaw'' (29)
LE = hfg(paKw paw'q ) (210)
where:
T7 : time averaged turbulent fluctuation of the air
temperature [K(]
u7 : time averaged turbulent fluctuations pf the horizontal
component of the wind velocity [ms"]
w : wind velocity in the vertical direction [m's1]
w : time averaged turbulent fluctuations of the vertical
component of the wind velocity [ms11
q : turbulent fluctuations in the specific humidity of the
air [kgvapor/kgdry air]
Assuming vertical gradients in temperature and absolute humidity
are insignificant as compared to the gradients in the direction as the
horizontal component of the wind, the latent and sensible heat can be
expressed as
H = CpaPaw'Ta (211)
LE = hfgPaw (212)
The equations are used as a basis for either the sensible or latent
heat flux in the atmosphere using measured values of the turbulent eddy
fluxes of heat and moisture. Specialized equipment was designed by the
Goddard and Pruitt to measure the parameters needed for the
calculations. However, the equipment was not reliable for low wind
speeds. This was attributed to the very small magnitudes of the
turbulent eddies for calm conditions.
13
Limitations of the mass balance (eddy flux and profile) methods
are the requirements for very sensitive equipment for measuring either
the vertical profiles or the vertical fluxes due to turbulence. King
(1966) questioned the practical use of the mass balance methods of
determination of the evaporation. Fritschen (1966) pointed out that
the mass balance methods do not integrate results over time, thus
requiring constant monitoring of meteorological parameters.
Another approach to estimate evaporation from the soil is the
energy balance method (Fritschen, 1966). The rate of change of
sensible heat in a control volume of air at the earth's surface is
equal to the rate at which sensible and latent heat are carried into
and out of the volume by the wind, the convection of sensible and
latent energy in the vertical direction, the energy transferred to the
soil and the net radiant exchange of energy. In the energy balance
method, it is assumed that the net transfer of sensible and latent heat
in the horizontal direction by the wind is negligible when compared to
the vertical movement of heat in the atmosphere. The equations of
motion are utilized as in the mass balance profile method with the
exception of using the similarity equation for the sensible heat flux
instead of the shear stress. A ratio of sensible to latent heat flux
(equations 24 and 22), referred to as the Bowen ratio (BR), is given
by
aT
H cpaPa Kh Jt
BR = (213)
8z
The energy balance at the soil surface is
Rn S G LE = 0 (214)
where:
Rn = net radiation incident upon soil surface [W*m2]
S = soil heat flux [Wm2]
G = sensible energy flux into atmosphere [Wm2]
LE = flux of latent energy from soil [Wm2]
An expression for the rate of evaporation of water from the soil can be
determined by substitution of equation (213) into (214) and
rearranging terms. This method is referred to as the Bowen ratio
equation and yields valid results for a wide range of conditions.
Fritschen states that it is imperative that efforts be undertaken to
assure that the assumption of no horizontal divergence of heat or
moisture be met in order for the Bowen ratio equation to yield
satisfactory estimates of the evaporation. Soil heat flux must be
measured as well, since omitting soil heat flux from the analysis could
lead to large errors.
In using the Bowen ratio, the soil surface temperature must be known as
in some of the other methods.
Penman (1948) used a combination of the energy and mass balance
methods with the objective of eliminating the need for surface
temperature. The underlying assumptions are the temperature of
evaporating surface is the same as the ambient air and the vapor
pressure is the saturated vapor pressure evaluated at the surface
15
temperature. It was also assumed that the net soil heat flux is zero
over a 24 h period. The resulting relationship
Pacpahfg (ea eas) + s(Rn G)
Ep = ( + s) (215)
where:
Ep potential evaporation [m3.m2]
ea = vapor pressure of atmosphere [kPa]
eas = saturated vapor pressure at the air temperature [kPa]
h = surface vapor transfer coefficient [m's1]
G = sensible heat flux in the atmosphere [Wm2]
Rn = net radiation upon the soil surface [Wm2]
s = slope of saturated vapor pressure line [kPa.'K1]
7 = psychrometric constant [kPaOK1]
provides an expression for the potential evaporation from the surface.
The boundary layer resistances account for advective conditions. The
Penman equation is the only method for estimating evaporation based
upon theoretical ideas and requires no highly specialized equipment.
Estimates of the evaporation have an accuracy of 5 to 10 percent on a
daily basis (Van Bavel, 1966). The disadvantage to the combination
method, as implemented by Penman, is that only daily estimates are able
to be determined and crop coefficients are necessary to estimate actual
evapotranspiration from various crops (Jones et al., 1984). The values
of the crop coefficients typically exhibit regional variation as well
as variation due to stage of crop growth.
Staple (1974) modified the Penman model to provide the upper
boundary condition for the isothermal diffusion of water in the soil.
The Penman model was modified by multiplying the saturated vapor
16
pressure by the ratio of vapor pressure at the soil surface to the
saturated vapor pressure at the same temperature (relative vapor
pressure). This incorporated the vapor pressure depression effect of
the soil water potential into the model. This seemed to match field
data for a clay loam soil. This model was not tested for more coarse
soils.
The Penman method eliminated the requirement of soil surface
temperature by neglecting the effects of soil heat flux at the expense
of being able only to predict the evaporation on a daily basis.
Various researchers have recently gone one more step in modeling the
evaporative loss of water from the soil by considering the movement of
water and heat below the soil surface (Van Bavel and Hillel, 1976;
Lascano and Van Bavel, 1983). Most of the research in which the
process of evaporation has been examined at this detailed level are
for a bare soil surface. This eliminated the complicating factors of
water removal by the plant, which surface temperature to use in
describing the driving potential for evaporation, and description of
the atmospheric boundary layer transfer coefficients.
Van Bavel and Hillel (1976) developed a model in which the partial
differential equations for water (liquid) and energy transfer in the
soil were considered to be independent processes. This approach
explicitly determined the soil heat flux for the boundary condition at
the soil surface as expressed in equation (214). The simultaneous
solution determined the actual evaporation from the soil as well as the
spatial distributions of water and temperature as a function of time.
The variation in hydraulic and thermal properties of the soil with time
was accounted for in this extended combination approach.
17
According to Fuchs and Tanner (1967) field and laboratory
observations indicated that the evaporation of water occurs in three
distinct stages. The first involves water evaporating from the soil
atmosphere interface. As the soil surface dries, the water must
evaporate at a location below the soil surface then the vapor diffuse
to the surface. The combination methods of simulating evaporation from
the soil reproduce this process fairly well. However, the evaporation
is assumed to occur at the soil surface. This has the effect of
lowering the soil surface temperature, when in the case of water
evaporating below the dry soil surface, the soil surface temperature
would be higher due to the reduced thermal conductivity of the soil.
This error in soil surface temperature was noted by Lascano and Van
Bavel (1983). Lascano and Van Bavel (1983 and 1986) verified the model
of Van Bavel and Hillel (1976) using soil water content and temperature
data obtained from field plots. During the earlier study, the
simulated soil surface and profile temperatures were found to agree
quite closely over the range of 25 to 37 oC while the model
underpredicted surface temperature by 2 to 5 OC when the soil
temperature was above 37 OC. Distribution of soil water is simulated
to within the expected error of measurement. Similar results were
obtained during the 1986 experiment.
Movement of water primarily in the vapor phase in the soil has
been observed by several researchers. Taylor and Cavazza (1954) noted
that the measured diffusion coefficients were higher than that for
water vapor in air and suggested that transport was due the combined
effects of convection and diffusion within the soil pore spaces.
Schieldge et al. (1982) used the following equations of the
18
conservation of mass and energy to simulate the diurnal variations of
soil water and temperature.
ae a a8 a aT aK
5 = z (DO ) + T (DT az) + T (216)
aT a aT
Cs at = ( )+ Q(O,T) (217)
where:
0 = volumetric water content [m3.m3]
t = time [s]
z = depth below soil surface [m]
T = soil temperature [OK]
Cs = volumetriq heat capacity of soil mixture
[Jemj3OKi]
S= thermal conductivity of soil mixture [Wm1OKI]
K = unsaturated hydraulic conductivity [m's1]
Do = hydraulic diffusivity [m2.s1]
Sao0 D a VS9Pw
S + PRT a80
Da = diffusivity of water vapor in air [m2.s1]
P = atmospheric pressure [kPa]
v = ratio of atmospheric pressure (P) to partial
pressure of air
S = soil porosity [m3m3]
g = acceleration due to gravity [ms1]
= soil matric potential [m]
Rw = ideal gas constant for water vapor
Pw = density of water vapor [kg'm3]
19
DT = diffusivity of water (liquid and vapor) due to
temperature gradients [m*'s10'1]
+ Da vSahon
Pw
Q(0,T) = heat flux within the spil due to movement of water,
liquid and vapor [W'm ]
a = soil tortuosity factor (= 2/3)
f= coefficient of thermal expansion [kg'm3"K1]
S= vapor transfer coefficient dimensionlesss]
h = relative humidity
These equations were originally developed and presented by de Vries
(1958) with good agreement with field observations being achieved.
Diurnal fluctuations of water content were damped out below a depth of
4 cm. Particular attention was given to the magnitudes of DT DTL and
DTV. It was noted that the DT was small during the night but has a
significant contribution to the flow of water during the day. The
vapor component (DT) was less than the liquid component (DTL) for
water contents greater than 0.30 m3/m3 and sensitive to changes in soil
temperature.
Jackson (1964) studied the nonisothermal movement of water using
equations (216) and (217) and concluded that classical diffusion
theory could be used satisfactorily if the relative vapor pressure is
greater than 0.97 with no modification to the diffusivity term.
However, if the relative vapor pressure is less than 0.97 then the
diffusivity must include the effect of the temperature gradient.
Model Objectives
The models presented thus far each have advantages and
disadvantages. The mechanistic models are based upon sufficient theory
that the diurnal fluctuations of the soil water and temperature
20
profiles can be observed. They also have the ability to provide
insight into the many processes involved in the evaporative loss of
water from the soil. Saxton (1986) stated that there is a need to
separate the evaporative loss of water from the soil from the losses
through the plant so that better understanding of the individual
processes can be achieved. The mechanistic models provide this
ability.
One might also note that in any of the models discussed
previously, conservation of water is discussed either in the vapor
phase, as in the mass balance methods, or the liquid phase, as in the
combination methods which consider the soil media. The vapor and the
liquid phases are not considered simultaneously.
The objectives in developing a new model were:
1. to account for the water changing from liquid to vapor
phase below the soil surface and diffusing to the
atmosphere;
2. to consider the overall mass continuity, specifically
include the liquid and vapor phases separately;
3. to account for the movement of water vapor in response
to temperature gradients in the soil; and
4. to simulate the diurnal variation of evaporation, soil
water and temperature profiles in response to conditions
at the soil surface.
Model Development
The soil for most purposes may be considered a continuous medium
in which the laws of conservation of mass and energy apply.
Application of the basic principles of thermodynamics to the soil
profile provide the basis for simulation of the temporal variation of
the distribution of water and temperature throughout the soil profile.
The following assumptions were made to simplify or clarify the
21
development of the model for coupled heat and mass transfer in the
soil.
1. The soil is unsaturated, therefore the soil water potential
is primarily due to osmotic and matric potential.
2. Water vapor behaves as an ideal gas.
3. The movement of water in the liquid and vapor phases occurs
due to concentration or potential gradients.
4. The liquid and vapor phases are in thermodynamic equilibrium
within the soil pore spaces.
Using the conservation of mass, the rate of change of the water
vapor within the soil air space was described in the model as
aPv a apv
t = (Dv  ) + Ev(z,t) (217)
where:
Pv = mass concentration of water vapor in the soil air
space [kgm3]
Dv = Diffusivity of water vapor in the soil air space
[m2.sl]
= Da
Da = Diffusivity of water vapor in the air [m2.s1]
a = soil tortuosity dimensionlesss]
Ev(z,t) = rate of phase change of water [kg'm3S1]
t = time [s]
z = depth below the soil surface [m]
The time rate of change of the concentration of water vapor at any
point in the soil space is a function of the rate of water vapor
diffusing from other regions of the soil and the rate at which water
changes from the liquid to the vapor phase. The evaporation rate is
22
considered to occur throughout the soil profile in response to
temperature and water content conditions in the soil. The rate of
phase change of the water (Ev) is considered positive if the water
changes from the liquid to vapor (evaporation), while condensation is
indicated by a negative phase change rate.
The diffusivity of water vapor in the soil air spaces is assumed
to equal that in the atmosphere, however, Schieldge et al. (1982), as
well as others, have made use of a tortuosity factor (a) to account
for the increased path length of the interstitial spaces of the soil
through which the water vapor must traverse. The tortuosity factor
included the effect of the local evaporation and condensation of water
vapor as it passes through the regions of differing temperature and
water content in the soil.
The volumetric water content is typically used to express the
amount of liquid water present in the soil. The mass concentration of
water in the liquid phase is obtained by multiplying the volumetric
water content by the density of water. The time rate of change of the
mass concentration of water is increased by a the rate of diffusion of
water from surrounding soil and decreased by the rate of evaporation.
The mathematical expression of the conservation of mass for the water
is
atW) = aQ EL(z,t) (218)
where:
0 = volumetric water content [m3m3]
Pw = density of water [kgm3]
QL = diffusion of liquid water [kgm2 s]
EL = rate of change of liquid water to water vapor
[kgmsils ]
For the range of temperatures generally occurring in the soil, the
density of water may be considered constant. The volumetric flow rate
of water (qL) can be obtained by dividing the mass flow rate (QL) by
the density of water (Pw). The volumetric flow rate of water in an
unsaturated soil is governed by Darcy's Law as follows:
qL = k z (219)
where:
qL = volumetric flow rate of liquid water [m3.m2s1]
k = unsaturated hydraulic conductivity [ms1]
S = soil water potential [m]
z = spatial dimension [m]
The relationship between soil water content and soil water potential is
unique for a given soil type. Therefore, the volumetric flux of water
can be expressed in terms of volumetric water content instead of soil
water potential by applying the chain rule of differentiation to
equation (219) as follows
L k a (220)
The hydraulic diffusivity of the soil is defined as
DL M K a (221)
a
By substituting the definition for the hydraulic diffusivity (221) and
the volumetric flow rate in terms of the soil water content (220) the
equation for the conservation of mass for liquid water (218) becomes
aB a a( E(z,t) (222)
aT = az (DL ) Pw
24
The conservation equations for the liquid (222) and vapor (217)
both have terms relating to the rate of change from the liquid phase to
the vapor phase. The vapor phase change ( Ev(z,t) ) is based upon the
volume of air while the liquid phase change term ( EL(Z,t) ) is based
upon the volume of soil. The phase change terms are related by the
fraction of soil volume occupied by the air which can be determined by
subtracting the volumetric water content from the total soil porosity.
The relationship between the vaporization terms is expressed as follows
E(z,t) = EL(,t) = (S O)Ev(z,t) (223)
Rearranging equation (223) and solving for Ev(z,t) yields
Ev(z,t) = Ez,t) (224)
Substitution of equation (224) into the vapor continuity equation
(217) results in the following equations describing the conservation
of water vapor within the soil.
SPv a aPv Ev(z,t)
= (D t ) (S ) (225)
The thermal energy equation describes the flow of heat within the soil
and includes the time rate of change of sensible heat in the soil
volume, the diffusion of heat due to temperature gradients, the
convection of heat due to diffusion of water and water vapor through
zones of variable temperature and the latent heat of vaporization. The
partial differential equation describing the rate of change of the
temperature of the soil used in the model was
aT 8 T 8 8apv aT
Cs " ( OF) + (PwCpwDL + cpvDvz ) zT
hfgE(z,t) j l (226)
where
X = thermal conductivity of soil mixture [WnmI.OI]
T = soil temperature [OK]
hfg = latent heat of vaporization [Jkg1]
Cs = volumetric heat capacity of soil mixture [Jm1K1]
= (1S)Pscps + Opwcpw + (SO)Pacpa
S = Soil porosity [m3m3]
ps = dry bulk density of soil [kg'm3]
Pa = density of moist air [kgm3]
cps = specific heat of solid portion of soil [JkgloK1]
Cpw = specific heat of water [Jkgl.'K1]
Cpa = specific heat of moist air [Jkg1oK1]
Equations (222), (225), and (226) describe the conservation of
energy and mass within the soil profile and account for the vapor and
liquid water phases separately. However, those three equations contain
the three state variables, soil temperature, volumetric water content,
and vapor concentration as well as the unknown rate at which water is
changed from the liquid to the vapor phase. Another independent
equation was needed to adequately describe the soil mass and energy
system.
If the water vapor were in equilibrium with liquid water, the
vapor would have a partial pressure corresponding to the saturated
vapor pressure and would be a function of the liquid water temperature.
Assuming that water vapor behaves as an ideal gas, then the saturated
26
vapor density can be determined from the saturated vapor pressure using
the ideal gas law.
pV = nRT
(227)
where:
p: pressure of the gas
V: volume of gas
n: number of moles of gas present
R: ideal gas constant
T: absolute temperature of gas.
The ideal gas law may also be applied for the individual components of
a gas mixture with the pressure being the partial pressure of the
component gas and using the number of moles of the component gas.
Substitution for the molecular weight of water and solving for the
density of water vapor one obtains
es(T)
R, T
where:
(228)
Pvs = concentration of water vapor t saturated vapor pressure
per unit volume of air [kgm" ]
es = saturated vapor pressure at given temperature [Pa]
T = temperature of free water surface [cK]
Rw = ideal gas constant for water vapor
= universal gas constant (R) divided by the molecular
weight of water
= 461.911 [m2.s2K1]
The above equation (228) is for the case in which the chemical
potential of the water is zero. However, due to matric and osmotic
27
forces, the water in the soil has a soil water potential less than
zero, thereby reducing the equilibrium vapor pressure from that of pure
water. The following relationship (Baver et al., 1972) can be used to
determine the saturated vapor pressure above a water surface with a
potential other than zero.
ev() = evs(T) exp( ) (229)
where:
ev(0) = saturated vapor pressure of water with chemical
potential, 0 [Pa]
= soil water potential [m]
g = acceleration due to gravity [mns2]
Using the ideal gas law, the vapor density over a water surface with a
chemical potential other than zero can be obtained by:
P = Pvs(T) exp( ) (230)
A set of simultaneous equations (222, 225, 226, 230) describe
the conservation of thermal energy, water vapor and liquid water for
the soil continuum and formed the basis of a coupled mass and energy
model for the soil. The assumption of thermodynamic equilibrium
between the liquid and vapor states yields the constitutive
relationship expressed in equation (230) and provides a fourth state
equation to be used in the model.
A wellposed problem also includes boundary conditions and, in
the case of transient problems, the initial conditions must also be
prescribed. The system of governing equations was onedimensional and
therefore required two boundary conditions. The first boundary is the
28
obvious soilatmosphere interface. It was assumed that the thermal
capacitance of the soil at the surface was negligible when compared to
the magnitudes of the fluxes which occur. Therefore, the net flux of
energy must be zero at the soil surface. The ability of the soil to
maintain a significant rate of evaporation at the soil surface was
assumed to be small as well. This assumption required the water to
evaporate at a finite distance below the soil surface rather than at
the soil surface. The boundary conditions for the energy, vapor and
liquid continuity were
aT I0 Pv
(PwcpwDL U + cpvDv  )Tav
"x  z=O
= Rn hh(Tz=0 Ta) (231a)
DL = P (231b)
apv
Dv = hm(pvO Pva ) (231c)
z=0
where:
cpw = specific heat of water [Jkg1.''1]
Cpv = specific heat of water vapor [Jkg1.K1]
DL = hydraulic diffusivity [m2s1]
Dv = diffusivity of water vapor in soil [m2s1]
hh = boundary layer heat transfer coefficient [Wm2.K1]
hm = boundary layer mass transfer coefficient [msl]
P = precipitation or irrigation rate [m3.m2s1]
Rn = net radiation incident upon soil surface [WNm2]
T = soil temperature [oK]
Tavg = average temperature [OK]
X = thermal conductivity of soil [W.ml1. 1]
Pva = mass oJ water vapor per unit volume of dry ambient air
[kgmJ]
PvO = concentration of water vapor at the soil surface, z=0
[kgm3]
0 = volume of water per unit total soil volume [m3.m3]
Ideally, for a semiinfinite medium, the flux of energy and mass
should be zero in the limit of depth (z) approaching infinity.
However, in anticipation of a numerical solution to the system of
partial differential equations, the boundary conditions were specified
at a depth of 1.5 meters. This depth was chosen by comparing the
error between the analytical and numerical solutions for conduction of
heat in a semiinfinite slab with constant uniform properties and a
uniform heat flux at the surface. This depth resulted in an error of
less than 0.1 OC at the lower boundary. The depth at which the
amplitude of the diurnal fluctuations in temperature is less than
0.1 OC is approximately 60 cm for a sandy soil (Baver et al., 1972).
The zero flux condition for the liquid and vapor continuity equations
represents an impermeable layer in the soil. This condition may or
may not physically exist in the field, but for most situations
encountered, the errors introduced into the solution at the depth of
chosen should be minimal. The boundary conditions used at the lower
boundary were
T = 0 (232a)
az Zzo
800
o I = 0 (232b)
a zz0Z
z = 0 (232c)
The system defining the movement of water vapor, liquid water and
energy is defined by a system of partial differential equations (222,
225, 226) with boundary conditions (231ac, 232ac). The energy
equation is coupled to the continuity equations by the rate at which
the water is changed from liquid to vapor and through the movement of
sensible heat associated with the flux of water between zones of
differing temperatures. The conservation of water vapor is coupled to
the soil temperature indirectly in the calculation of saturated vapor
pressure. This coupling, along with the variation of the soil
properties with time and space, renders an analytical solution beyond
reach, therefore; a numerical solution was required.
Determination of Model Parameters
Solution of the governing equations for the energy and mass
balance for the soil requires knowledge of the properties relating to
the soils ability to diffuse heat, water, and water vapor. The other
parameters to be determined relate to the rate at which heat and water
vapor are dissipated from the soil surface to the air. In order to
determine a solution, relationships for determining thermal, hydraulic
and vapor diffusivities as well as the surface transfer coefficients
must be determined.
Diffusivity of Water Vapor
Diffusivity is the constant of proportionality relating flow to
the gradient in potential as stated in Fick's Law of molecular
diffusion. For the case of water vapor diffusing through air, the
31
diffusivity, or coefficient of diffusion, is a measure of the number of
collisions and molecular interaction between molecules of water vapor
and the other constituent components of the air (ASHRAE, 1979) and is a
function of temperature, total air pressure and partial pressure of
water vapor. Eckert and Drake (1972) presented the following equation
to calculate the diffusivity of water vapor in air.
Dva 2.302 ( (T X 105 (233)
va p
where
Dva = diffusivity of water vapor in air [m2.s1]
p = atmospheric pressure [Pa]
Po = reference atmospheric pressure [Pa]
= 0.98 X 105 Pa
T = air temperature [OK]
T = reference air temperature [OK]
= 256 OK
However, for water vapor in the soil air space, the path is more
convoluted, thus increasing the probability of collisions with other
particles and lowering the kinetic energy of the water vapor molecules.
To account for the increased path length within the soil, a tortuosity
factor (a) has been introduced to reduce the effective coefficient of
diffusion in other porous materials. De Vries (1958) used a value for
the tortuosity of 0.667. Using this information the coefficient of
water vapor in the soil (Dv) becomes
Dv = Dva
= 1.535 (P T 1.81 X 105 (234)
~p ) To)
For the purposes of simulation, the barometric pressure (p) was assumed
to equal the standard pressure (po). The diffusivity then became a
function of time and space due to the temporal and spatial variation of
the soil temperature. Diffusivity was determined throughout the soil
profile by substitution of the soil temperature in equation (234).
Thermal properties
Very little literature exists regarding measurement of the thermal
properties of the soil. Soil composition as well as density and water
content affects the thermal properties of the soil (Baver, 1972). The
thermal diffusivity is the ratio of thermal conductivity to the product
of soil density and specific heat. Density is a property which can be
obtained as a function of depth at a given location by core samples.
Vries (1975) describes a method by which the volumetric heat
capacitance and the thermal conductivity of the soil can be calculated
based upon the volume fractions of the various soil constituents.
The volumetric heat capacitance is defined as the product of the soil
density and the specific heat and can be calculated from:
Cs = qCq + XmCm + XoC + ww + xaCa (235)
where
C = volumetric heat capacity [Jm3.K1]
x = volume fraction [m3.m3]
q = quartz
m = mineral
o = organic
33
w = water (same as volumetric water content)
a = air
s = soil composite
The difficulty arises in estimating the various volume fractions of the
different components. For quartzitypic sands, approximately 70 to
80 percent of the solid constituent of the soil is quartz, with
generally less than 1 to 2 percent organic material and the remainder
consisting of other ninerals. According to information presented by
DeVries (1975), the volumetric heat capacity and density of quartz and
other minerals are very similar. For the purposes of this simulation,
the volume fraction of the organic material was assumed to be zero.
The thermal conductivity cannot be calculated in such a straight
forward manner. Thermal conductivity is defined as the constant of
proportionality in relating the heat flux by conduction to the
temperature gradient. Conduction of heat occurs due to physical
contact between adjacent substances. In a solid material such as steel
or concrete which is fairly homogeneous, the material conducting heat
can be considered continuous. However, the soil is a mixture of solid,
liquid and gaseous components and the area of physical contact between
soil particles is a function of the particle geometry. When the soil
is dry, the area of contact may be a single point. As the soil is
wetted, a thin film of water adheres to the soil particle and increases
the contact area between adjacent particles. This increase in area
accounts for a rapid change in the thermal conductivity as the soil is
initially wetted (Figure 21). Over a range of water contents between
dryness and saturation, the increase in thermal conductivity becomes
linear with increasing water content. As the soil approaches
saturation, the thermal conductivity begins leveling off to some
relatively constant value. Another complicating issue in calculating
the thermal conductivity of the soil is that the temperature gradient
in the solid and liquid portion of the soil can be significantly
different from that in the gaseous phase of the soil.
De Vries (1975) presented a method by which the thermal
conductivity could be calculated using a weighted average of the
various soil constituents where the weighting factors were a product of
the individual volume fractions and geometric factors as follows
kqXq + km xmA + k oX 0 + k X + kaxaa (236)
Sqxq + kmXm + koXo + kX + kaxa
where
S = thermal conductivity of the soil [Wm1K1]
\ = thermal conductivity of soil component [W'm1.KI1]
xi = volume fraction of soil component [m3m3]
ki = geometric weighting factor dimensionlesss]
i = quartz (q), mineral (m), organic (o), water (w), air (a)
The geometric weighting factor depends upon geometric
configuration of the soil particles and the incorporated void space.
It represents the ratio of the spatially averaged temperature gradient
in the ith soil component to the spatially averaged temperature
gradient of the continuous component of the soil (usually water). For
example, kq represents the space average of the temperature gradient
in the quartz particles in the soil to the space average of the
temperature gradient in the water. The geometric weighting factor can
be determined by
ki (kia + kib + kic) (237)
where
kij = the ratio of the temperature gradient in the ith
component to the temperature gradient in the continuous
component in the direction of the jth principal axis
i = quartz (q), mineral (m), organic (o), water (w), air (a)
j = a, b, or c for each of the three principal axis of the
particle
The ratio of the temperature gradients in the ith component can be
determined by the following:
kj = 1 + 1 (238)
+ ( x 1
The shape factor for each principal axis (gj) can be approximated by
various empirical relationships depending upon the ratio of the unit
vectors (ua, ub, uc) of the principal axes of the soil components
(Table 21). The sum of the three shape factors must be unity.
In most cases, water is considered to be the continuous phase of
the soil in determining soil thermal conductivity. However, as the
soil dries and the film adhering to the surface of the soil particle
begins to break, making air the continuous phase. Equation (237) can
be used in these cases replacing the thermal conductivity of water (Aw)
with the thermal conductivity of air (Xa). De Vries (1975) noted that
the values for the thermal conductivity in the case of air being the
continuous phase were consistently low by a factor of approximately
36
1.25. This method yielded values of thermal diffusivity within ten
percent of those measured (de Vries, 1975). Extensive detail regarding
calculation of the thermal conductivity of the soil is given in
de Vries (1963).
The thermal properties were calculated using equation (235) for
the volumetric heat capacity and equations (236), (237) and (238) to
determine the thermal conductivity. Volume fractions of the various
soil constituents was determined from soil classification data and
knowledge of the soil bulk density and porosity. Shape factors used
for calculation of the thermal conductivity were for a typical sand
grain (Table 21).
Hydraulic Properties
The parameter governing the movement of liquid water in the soil
is a measurable property of the soil and is analogous to the thermal
diffusivity. The hydraulic diffusivity is a derived property of the
soil (i.e. not directly measured) and is defined as the hydraulic
conductivity divided by the specific water capacity of the soil (Baver
et al., 1972). The hydraulic conductivity is the constant of
proportionality for the diffusion of water in response to a gradient in
the soil water potential (Figure 22), while the specific water
capacity is the slope of the soil water retention curve. The
hydraulic conductivity (Figure 23) and the specific water capacity
(Figure 22) vary according to the soil composition as well as the soil
water potential. The measurement of hydraulic conductivity can be
accomplished by several methods, but most all require meticulous
control of the potential gradients and a great deal of time. This is
especially true if measurements are desired over a wide range of soil
37
water potential. Measurements of the soil water retention curve
require considerably less detail and are generally published for a
wide variety of soil types. Van Genuschten (1980) presented a method
by which the hydraulic conductivity could be calculated for unsaturated
soils using the soil water retention curve and the saturated hydraulic
conductivity. The Van Genuschten approach involves estimation of an
equation for the soil water retention curve of the form
9 = ( 10 r 1) m (239)
Ts r 1 + (a)n
where
9 = dimensionless water content
0 = volumetric water content [m3m3]
r = residual volumetric water content [m3m3]
es = volumetric water content at saturation [m3m3]
= soil water potential [ m ]
m,a = regression coefficients
n = (1 m)
The coefficients m and a are determined by nonlinear least squares
regression of the soil water retention curve. That function is then
substituted into equations for the hydraulic conductivity presented by
Mualem (1976). The resulting expression for the diffusivity is
D(B) = (1 m) k (0.5 I/m) [(1e1/m)m+ (1_01/m)m_2] (240)
am(OsOr)
The Van Genuschten method yields a continuous function for the
hydraulic diffusivity which is highly desirable for numerical
38
simulations over the range of the water contents expected to occur in
the fields.
For modeling purposes, the Van Genuschten method was employed for
published potentialwater content data for a Millhopper fine sand
(Carlisle, 1985). The water release curve usually varies with depth
due to the spatial variation of soil composition. Data for the a
single water release curve for a uniform soil profile was obtained by
averaging the volumetric water content for the specified water
potential over the Ai and A2 horizons of a Millhopper fine sand. The
nonlinear regressions were then performed to yield the residual water
content, and the regression coefficients, a and n.
Surface Transfer Coefficients
The final parameters to be estimated are the surface heat and mass
transfer coefficients. Penman (1948) used an empirical wind function
to calculate the mass transfer coefficient based upon wind speed.
hm' a WSb (241)
where
hm = surface mass transfer coefficient
WS = wind speed
a,b = empirical constants
Another approach is to use the equations of motion to describe mass and
energy transfer within the atmospheric boundary layer. Brutsaert
(1982) presented a detailed review of the equations of motion as
applied to the atmosphere. The basic assumption in these analyses is
that the boundary layers for momentum, energy and mass are similar.
Consequently, the surface sublayer becomes the area of most concern and
39
is defined as the fully turbulent region where the vertical turbulent
fluxes of mass and energy do not change appreciably from that at the
surface (Brutsaert, 1982). In other words, the vertical flux of mass
and energy is constant.
According to Brutsaert (1982), Prandtl introduced the use of the
logarithmic wind profile law into meteorology in 1932. This is an
approximation of the wind velocity profile in the surface sublayer
WS = In (0) (242)
where
WS = wind speed [ ms1]
U* = shear velocity
1/2
= TI
P
z0 = surface roughness height [ m ]
k = von Karman constant dimensionlesss]
T = shear stress at the surface [Nm21
p = density of air [kgm3]
According to Sutton (1953), for practical purposes, the shear velocity
can be estimated as
WS
U = (243)
The bulk mass transfer within the surface sublayer is defined by
Qm = CeyWSr (Pvs Pvr) (244)
where
Qm = vertical mass flux [kg'm2.s1]
WSr = wind speed at reference height, zr [m's1]
Cer = Dalton Number; dimensionless mass transfer coefficient
Pvr = water yapor concentration at reference height, zr
[kg'm ]
Pvs = water vapor concentration at soil surface [kg'm3]
The surface mass transfer coefficient (hm) used in this model is
related to Cer by
hm = WS, Cer (245)
Brutsaert (1982) presented functions for the Dalton number in terms of
the drag coefficient and the dimensionless Schmidt number as
1/2
Cd,
Cer = (246)
1/2
(B + Cdr)
where
B = function of dimensionless Schmidt number
Cer = Dalton number
Cdr = surface drag coefficient dimensionlesss]
2
WS
41
The empirical function, B, depends upon whether or not the surface is
hydrodynamically smooth or rough. In general, if the Reynolds number
based upon the shear velocity and the roughness height is less than
0.13 then the surface is considered smooth and B is determined using
Equation (247). The surface is rough if the Reynolds number is
greater than 2.0 and Equation (248) is used to determine the value
of B.
2
/3
B = 13.6 Sc 13.5 (247)
/4 2/
B = 7.3 Re 4 Sc 5.0 (248)
where
Sc = Schmidt number
Dv
Re = Reynolds number
z0 U*
The dimensionless heat transfer coefficient or Stanton number
(Chr) can be determined by substitution of the Prandtl number for the
Stanton number for the Schmidt number in equations (247) and (248)
above. The surface heat transfer coefficient can then be determined by
hh = Chr Pacpa WSr
(249)
42
Fuchs et al. (1969) used a method in which instability in the
lower atmosphere was accounted for in conjunction with the logarithmic
wind profile to calculate the surface mass transfer coefficient. This
involved using the KEYPS function (Panofsky, 1963) to determine the
curvature of the diabetic lapse rate (w) of the logarithmic wind
profile. The surface mass transfer coefficient was calculated from
2
h = k2WS ( 7 + In (z + D) (250)
m Zo
where
hm = transfer coefficient of mass from surface to height, z,
in the air [ms1]
k = von Karman constant dimensionlesss]
WS = wind speed at height, z [ms1]
D = height displacement [m]
= d + z0
d = height above soil surface where wind velocity is zero
[m]
Z0 = roughness length [m]
z = height above the displacement height [m]
x = curvature of the diabetic lapse rate or diabetic
influence function dimensionlesss]
Fuchs et al. (1969) then calculated the surface heat transfer
coefficient using
hh = hmcpaPa
(251)
The diabetic influence function accounted for the transfer of air
movement in the vertical direction due to density gradients caused by
temperature gradients and is a function of the Richardson number
(Figure 24). Fuchs et al. (1969) stated that the roughness height
varied from 0.2 to 0.4 mm for a bare soil surface and accounted for a
small variation in the calculated heat transfer coefficient.
Therefore, for the purposes of this study, a roughness length of 0.3 mm
will be used Fuchs et al. (1969) also noted that the height of zero
wind velocity (d) was zero for a bare soil surface. This term was
employed in the approximate wind profiles to account for the fact that
wind does not penetrate full vegetative canopies and for practical
purposes the surface where the wind velocity is zero occurs at some
finite height above the soil surface (Brutsaert, 1982; Sutton, 1953).
Fuchs et al. (1969) compared transfer coefficients calculated using
equation (250) to that determined from field data for a bare sandy
soil and obtained fairly close agreement.
For the purposes of this model, the approach used by Fuchs et al.
(1969) to determine the surface mass (Equation 250) and heat
(Equation 251) transfer coefficients was employed. The diabetic
influence function was utilized to account for atmospheric instability
as proposed by Fuchs et al. (1969). A roughness height (Zo) of 0.3 mm
was utilized. The equations used in the formulation of this model and
the determination of parameters are summarized in Table 22
Numerical solution
The partial differential equations used to describe the mass and
energy balance in the soil must be solved numerically since analytical
solutions are not possible for the coupled nonlinear equations. Many
44
solution techniques utilize the finite difference form of the
differential equations. The differential operator in the continuous
domain becomes a difference operator in the discretized domain.
Difference operators are either forward, backward, or central
difference operators. The difference operators for the function of the
independent variable, x, [ f(x) ] are
af f(xj+l) f(xj)
forward: af = Af f(x ) f(xj) (252)
ax xj+1 xj
af f(xj) f(xj1)
backward: = Vf = (253)
Xj xj1
af
central: = 6f (254)
1 f(xj+1) f(xj) f(xj) f(xj1)
T xj+1 xj xj xj1
In expressing the differential equations in difference form, the
derivative with respect to time is accomplished using the forward
difference operator, while spatial derivatives are expressed using any
of the three difference operators. The continuous domain must first be
discretized or divided into several discrete regions such that the
finite difference approximation of the differential equation approaches
the differential equation in the limit of the grid spacing going toward
zero (Figure 25).
The partial differential equation defining the conservation of
energy within the soil profile (equation 226) stated that the change
in sensible heat over time is due to heat transferred by conduction
plus sensible heat carried by the diffusion of water in the liquid and
45
vapor phases from a region of differing temperature less the latent
energy required for evaporation of water. Examining this for the jth
node in the discretized domain, the above equation can be expressed in
terms of heat flux transferred (q) into node j from the surrounding
nodes as.
dxdydwj)qsj = dxdy (qcj+l,n + qcjl,n + qLj+l,n
+ qLjl,n + qvj+l,n + qvjl,n)
(dxdydwj)qej,n
(255)
The subscripts
s :
c :
e
L :
v
j :
j+l
j1
n
in the above representation denote
change in sensible heat
heat conduction
latent heat
heat transfer due to diffusion of water (liquid)
heat transfer due to diffusion of water vapor
the current node
the node immediately following node j, and
the node immediately preceding node j.
the current time step
and
dxdy :
dxdyddwj :
dwj :
The change in
to n+1 is
crosssection area normal to z axis
volume of the node j
height of the cell for node j
sensible heat per unit volume of node j from time step n
Tj,rn+ Tj,
qsj = Cs dt
(256)
The heat flux conducted from nodes j+1 and j1 to node j at a given
time step, n, is described by Fourier's law of heat conduction and is
expressed as follows
Xj
qcj+l,n= (Tj+i,n Tj,n) (257)
qcjl,n= (Tjl,n Tj,n) (258)
dzj1
noting that dzj represents the difference between the depth of nodes
j+1 and j.
In developing the differential equations for this model, it is
assumed that mass movement is due to diffusion and obeys Fick's law of
diffusion. Therefore, the heat transferred by diffusion of water and
water vapor from node j+1 to j is expressed as
S(Oj+,n Oj,n) (Tj+1,n+ Tj,n)2
qLj+l,n = PwCpwjDLj dz2 (259)
(Pvj+,n Pvj,n) (Tj+1,n+Tj,n)
qvj+1,n = cpvjDvj dz 2 (260)
and similarly for the heat flux transferred by mass diffusion from node
j1 to node j as
(0jl,n 1j,n) (Tjl,n+Tj,n)
qLjl,n = PwjCpwjDLj dzj1 2 (261)
S(Pvjl,n Pvj,n) (Tjl,n+Tj,n)
qvjl,n = cpvjDvj dzj 2 (262)
The latent heat associated with the phase change of water for the
current time step (Ej,n) must be transferred per unit volume of soil
associated with node j and is determined by
hfgjEj,n
(263)
qej,n =
Substitution of the numerical expressions for the individual terms of
the energy balance (equation 255) and dividing by the crosssectional
area normal to the z axis (dxdy) yields
Tj,n+lTj,n Tj+l,nTj,n Tj1,n Tj,n
dwjCs dt= J dz + dzj_
(Oj+1,n0j,n) (Tj+l,n+Tj,n)
+ PwjCpwjDLj dzj 2
(0ji,n 0j,n) (Tjl,n+Tj,n)
+ PwjcpwjDLj dzj 2
(Pvj+1,n Pvj,n) (Tj+1,n+Tj,n)
+ cpvjDvj dzj 2
(Pvjl,n Pvj,n) (Tjl,n+Tj,n)
+ cpvjDvj dzj. 2
dwj (hfg,jEj,n) (264)
The partial differential equation describing the diffusion of
water within the soil can be transformed to difference form in a
similar fashion as the energy equation. The change in water content
of a volume of soil for the jth node is caused by diffusion of water
from nodes j+1 and j1 to node j less the amount of water changed to a
vapor phase. The numerical expression for the water continuity per
unit area becomes
48
(0j,n+l Oj,n) 0 +1,n 0j,n (Ojl,n Oj,n)
dwJ dt DLj dzj + dzjl
Ej,n
dwj (265)
Pw,j
Similarly, the differential equation describing the conservation of
mass in the vapor state is transformed as
(Pvj,n+1 Pvj,n) (Pvj+l,n Pvj,n)
S dt Dvj dzj
(266)
(Pvjl,n Pvj,n) dwj Ej,n
+ Dv dzj1 + (Sj j,n)
Equations (264), (265), and (266) constitute the numerical
equivalents for the partial differential equations describing the
conservation of thermal energy and mass in the soil profile. However,
the conservation relationships have yet to be developed for the
boundaries at the soil surface and the bottom of soil profile. Since
there is no volume associated with the surface node (Fig. 21), there
is no storage capacity for energy or mass. Therefore, the flux of mass
or energy from the node directly below the surface (j=2) to the surface
added to the net flux due to solar radiation onto the surface and the
heat flux due to convection from the air to the surface must be zero.
This can be expressed directly in difference form as follows.
energy: (T2,n T,n) = Rnn + hhn(Tan T,n)
energy:  (T2,n T, ) = Rn, + hhn(Tan T1,n)
+ PwlcpwlDL1
DL1
water: 
Dvi
vapor: T
(267)
(02,n 1l,n) (Pv2,n Pvl,n) (T2,n + Tl,n)
dzl + CpvlDv1 dz1 2
(02,n 91,n) + P = 0
(Pv2,n Pvl,n) + hm(Pva Pvl,n) = 0
(268)
(269)
: net solar radiation incident
at time step n [Wm'2]
: precipitation at time step n
: boundary layer mass transfer
step n [ms' ]
upon the soil surface
[m3.m3.s1]
coefficient at time
: boundary layer heat transfer coefficient at time
step n [Wm 2. ]
: air temperature at time step n [JK]
: ambient water vapor concentration at time step n
[kgm3]
The boundary conditions at the lower boundary of the soil are
developed in a similar manner. The last node has volume and
capacitance for storage of energy and mass. A zero flux boundary
condition is used. There are two ways in which the zeroflux condition
can be represented. These will be developed for the water for the
purpose of illustration.
where:
P
hmn
hhn
Tan
Pvan
50
If the water content of the last node (j=nc) were equal to that of
the preceding node (j=nc1), then no water would flow between the two
nodes. This would cause both of the last two nodes to act as one node
during the simulation. The other viewpoint would be that the sum of
the flows into the last node from the preceding node and an imaginary
node following must be zero.
DLnc DLnc
Sz (On+c+l,n nc,n) =  l(ncl,n Onc,n) (270)
If the distance between nodes nc and nc1 is the same as the distance
between the nodes nc and the imaginary node nc+1, then it follows that
Onc+l,n = Onc1,n (271)
Substitution of 271 into 265 yields as the boundary condition for the
node, j=nc, the following
(Onc,n+1 Onc,n) DLnc
dwnc dt 2 d (Onc1,nnc,n)
(272)
nc,n
dwnc wnc
Pwnc
The same approach produces the following as the boundary condition for
the vapor continuity and energy equations, respectively.
(Pvnc,n+l Pvnc,n) Dvnc
dwnc dt = 2 (Pvnc1Pvnc)
(273)
dwnc Enc,n
(S Oj,n)
(Tnc,n+1Tnc,n) (Tncl,n Tnc,n)
dwnc Cs,nc dt = 2 fnc dznc1
(nc1,n Onc,n) (Tncl,n+ Tnc,n)
+ 2 PwncCpwncDLnc dznc1
(Pvncl,n'Pvnc,n) (Tncl,n+ Tnc,n)
+ 2 CpvncOvnc dznc 2
+ dwnc Enc,n (274)
The conservation relationships provide sufficient information to
determine three of the desired quantities for the soil profile leaving
a fourth remaining unknown. The constitutive relationship requiring
that the water in the liquid phase be in thermodynamic equilibrium with
the vapor phase was used to provide the remaining equation. Since the
surface node provides only an interface between the soil and the
atmosphere, this equilibrium condition was not necessary for the
surface node. Under atmospheric conditions, the soil may not become
completely dry (0 % volumetric water content), but will reach some
moisture content which is in equilibrium with the atmosphere. This is
typically taken to be the same as the permanent wilting point of the
soil (0= 150 m) and is the same as the residual water content
described by van Genuschten (1980). As the soil surface dries, the
soil reaches this equilibrium moisture content and is assumed to act as
an interface between the air and the wet soil below.
The system of equations (264) through (269)and (272) through
(274) along with equation (230) represent an explicit solution to the
differential equations, that is, the values of the state variables at
the next time step are functions of those at the previous time step.
Using the explicit representation allows a simple algorithm to be used
52
during the solution phase of the system of equations and is equivalent
to an Euler integration in time (Conte and de Boor, 1980). The
disadvantage of the explicit solution is that the numerical error may
propagate through time and grow. The maximum time step is functionally
related to the boundary conditions and thermal properties of the soil.
The system of equations can also be changed such that all
derivatives are evaluated at the next time step in which none of the
state variables are known. This implicit representation requires a
relatively complicated iterative or matrix solution technique for the
system of equations. The advantage of implicit solution methods is
that the value of the state variables do not depend upon previous
values implying that the only source of error would be due to roundoff
or truncation errors and would not propagate or grow with time and
would yield an unconditionally stable solution (Conte and de Boor,
1980).
Another technique which has some of the desirable characteristics
of both the explicit and implicit methods is the alternating direction
(ADI) method. This is accomplished by incrementally marching through
space in one direction (z=0 to zo) evaluating the derivatives
containing the previous node (j1) at time, t + 1/2 dt, and those
containing the following node (j+1) at time, t. Then returning in the
opposite direction (z=zo to 0), evaluate the derivatives containing the
node j+1, at t=t+dt, and the derivatives containing j1 at t=t+1/2dt.
This technique uses a relatively simple algorithm similar to that for
the explicit methods because all of the values of the state variables
used in estimating the state variable at the next time step are known
quantities. However, the number of equations to be evaluated is twice
53
that of the standard explicit methods. Increased stability is obtained
over the standard explicit methods but is less than that for the
implicit techniques. This implies that using the same grid spacing, a
larger time step can be used in the ADI methods than that for explicit
methods.
An ADI finite difference technique was used to solve the system of
equations for the soil profile due to the increased stability
characteristics over explicit methods. Appendix A contains a detailed
development of the numerical equations used in the ADI technique.
Hydraulic (equation 240) and thermal properties (equations 235 and
236) of the soil and the surface transfer coefficients (equations 250
and 251) were calculated at the beginning of each time step. A
general description of the solution algorithm is shown in Figure (26).
The energy and water equations were solved for the temperature and
volumetric water content, respectively at the next time step. The
equilibrium condition was used to determine the vapor concentration by
substitution of the values of soil temperature, water content, water
potential, and saturated vapor concentration in equation (230). The
vapor continuity equation was used to determine the evaporation rate.
The calculations were repeated using the new evaporation rate
until the absolute value of the maximum fractional change in any of the
variables was less than a prescribed convergence criterion.
The system of numerical equations was solved using computer code
written in FORTRAN77. A variable grid spacing was used throughout the
soil profile with the smaller mesh being located near the soil surface
due to expected large gradients in temperature and moisture content
once the soil surface begins to dry (Table 23). The program was then
54
run to simulate the response of the soil system under constant
radiation, ambient air temperature and relative humidity conditions for
a period of 72 hours using various time steps to evaluate the numerical
stability and numerical error characteristics of the procedure.
Thermal and hydraulic properties of the soil were assumed to be
constant for the duration of the simulation for an Millhopper fine
sand. The water content and temperature profiles were assumed to have
a uniform initial distribution as were density and soil porosity.
Overall mass and energy balances were calculated for the soil
profile during the simulation. Any residual in these balances would
constitute error arising from the numerical algorithm and was
accumulated on an hourly basis for analysis purposes. The sum of the
squares of the residuals for the energy and mass balances for the
simulation period are presented in Figures 27 and 28, respectively,
for the different time steps. The sum of the squares of the residuals
(SSR) for the energy balance decreased when the time step of was
increased from 30 to 60 s then remained fairly constant for larger time
steps. The square root of the SSR represents an estimate of the
standard error of the estimate of the total heat flux. For the time
steps of 60, 120, 300 and 600 s, this represents approximately one
percent of the cumulative soil heat flux. The standard error of the
mass balance was approximately two percent of the cumulative
evaporation over the 72 h simulation (600 s time step). It is expected
that as properties are varied with time, that the standard error would
increase due to increased gradients in soil temperature and water
content particularly near the surface. Therefore, a time step of 60 s
was utilized in the model validation and sensitivity analysis.
Table 21.
Object
Formulae for the geometric shape factors used in
calculating thermal conductivity based upon ratios
of the unit vectors of the principal axes.
Sphere
Ellipsoid of revolution
x=y=nz n = 0.1
n = 0.5
n = 1.0
n = 5.0
n = 10.
Elongated cylinder with
elliptical xsection
x = ny;
0.14 0.14
0.33
0.49
0.41
0.33
0.13
0.07
(n + 1)1
0.33
0.49
0.41
0.33
0.13
0.07
n/((n +1)
0.33
0.02
0.18
0.33
0.74
0.86
0.00
0.72
Typical sand grain
Table 22. Summary of equations used in a model of heat and
mass transfer in the soil to determine state
variables and transport parameters.
Process or Parameter Equation No.
Conservation of mass (liquid) 222
Conservation of mass (vapor) 225
Conservation of energy 226
VaporLiquid equilibrium 230
Surface Boundary Conditions
Conservation of mass (liquid) 231b
Conservation of mass (vapor) 231c
Conservation of energy 231a
Lower Boundary Conditions
Conservation of mass (liquid) 232b
Conservation of mass (vapor) 231c
Conservation of energy 232a
Soil properties
Diffusivity of water vapor in soil (Dv) 234
Volumetric heat capacity (Cs) 235
Thermal Conductivity (X) 236
Hydraulic Diffusivity (DL) 240
Surface Mass Transfer Coefficient (hm) 250
Surface Heat Transfer Coefficient (hh) 251
Table 23.
Mesh spacing used in grid generation for the coupled
heat and mass transfer model for sandy soils.
Cell Width
(cm)
0.0
1.0
1.0
1.0
1.0
1.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
2.0
5.0
5.0
5.0
5.0
5.0
5.0
5.0
10.0
10.0
10.0
10.0
10.0
10.0
10.0
Distance to
next node
(cm)
0.5
1.0
1.0
1.0
1.0
1.5
2.0
2.0
2.0
2.0
2.0
2.0
2.0
3.5
5.0
5.0
5.0
5.0
5.0
5.0
7.5
10.0
10.0
10.0
10.0
10.0
10.0
Node
Depth
(cm)
0.0
0.5
1.5
2.5
3.5
4.5
6.0
8.0
10.0
12.0
14.0
16.0
18.0
20.0
23.5
28.5
33.5
38.5
43.5
48.5
53.5
61.0
71.0
81.0
91.0
181.0
191.0
201.0
58
15
12
E
E
0
o 3
0)
0.0 0.1 0.2 0.3 0.4 0.5 0.6
Volumetric Water Content (m3.m3)
Figure 21. Thermal conductivity of a typical sandy soil as a function
of volumetric water content (DeVries, 1975).
59
104 0.010
Potential
1 \ _ Capacity / 0.008 E
E 103 
/ ooo8
S102
0
0.004 )
o /
10
0 /
L I 0.002 L
/ 0
1 0.000
1 I I I 0.000
0.00 0.10 0.20 0.30 0.40
Volumetric Water Content ( cm3 / cm3 )
Figure 22. Soil water potential (tension) and specific water capacity
as a function of volumetric water content for a sandy soil.
60
102
10  Diffusivity 10
~1 Conductivity o105
E 104
1 0 103
102 
0 102
o 104..
105101
o 10
106 /*102 1
1 10 ..103
108 I 10o4
0.00 0.10 0.20 0.30 0.40
Volumetric Water Content ( cm3/cm3)
Figure 23. Hydraulic conductivity and hydraulic diffusivity as a
function of soil water content for a typical sandy soil.
1.50 1.00 0.50 0.00
0.50
0.50
Richardson Number ( Ri)
Figure 24. Relationship of the diabetic influence function and the
Richardson number.
0.5
0.01
1.5
2
2.00
I i ,
0.5
1.0 
H=f(To,hh,WS)
P Rn Qv=f(To,RH,hm,WS)
j=2
j=3
AZ j jl
i+l
ncI
nc
Figure 25. Schematic of discretized domain for evaporation and soil
temperature model.
Figure 26. Flowchart describing the solution algorithm to solve the
system of finite difference equations for a coupled heat
and mass transfer model.
3.
0
2
o
S2
E
4,
0
o
0
o" 1
E
C
0
60 120
Time Step ( s )
300
600
Figure 27. Sum of squared residuals for the soil profile energy
balance after a 72 hour simulation using various time
steps.
1
_1_1111 _
Sq
Sq
60 120
Time Step (
K
300
600
Figure 28. Sum of squared residuals for the soil profile mass balance
after a 72 hour simulation using various time steps.
0.06
0.04
0.02
if fll
m\c~
i
T
CHAPTER III
EQUIPMENT AND PROCEDURES FOR MEASURING THE MASS AND
ENERGY TRANSFER PROCESSES IN THE SOIL
Introduction
The process of evaporation of water from the soil is a complicated
process in which the dynamic processes of heat and mass transfer are
interrelated. Lysimetry has been defined as the observation of the
overall water balance of an enclosed volume of soil. Lysimeters have
been used for over 300 years to measure soil evaporation and crop water
use (Aboukhaled et al., 1982) and can be generally classified as either
weighing or nonweighing. As their names imply, weighing lysimeters
measure evaporation by monitoring the weight changes within the
enclosed soil volume while the nonweighing, or drainage lysimeters
determine evaporation by monitoring amount of water drained from the
bottom of the tank and the amount of rainfall upon the lysimeter. Most
drainage lysimeters are monitored on a seven to ten day cycle while the
weighing lysimeters can be monitored continuously to obtain hourly
evaporation rates.
The energy status of the soil can generally be determined by
measuring the temperature distribution over time. The energy flux at
the surface must be determined as well as latent heat transfer to
facilitate a complete energy analysis. Surface heat flux occurs as
radiant and convective heat transfer. The radiant energy impinging
upon the soil surface is either absorbed or reflected. The amount of
reflected and absorbed energy is primarily dependent upon the soil
67
surface. The reflectance and absorbtance change dramatically as a
function of water content. Convective heat losses from the surface are
generally estimated rather than measured using empirical functions to
determine surface heat transfer coefficients.
Objectives
The processes of energy and mass transfer in the soil are most
readily represented by the rate at which water is lost from the soil
(evaporation rate), the cumulative water loss (cumulative
evaporation), temperature, volumetric water content, and water vapor
density. All are temporal in nature, while the latter three are
functions of depth as well. The goal of the experimental design was to
provide data for the validation and calibration of the transport model
described in the previous chapter. Specific objectives were
1. to design a lysimeter and calibrate the instrumentation to
measure hourly evaporation rates,
2. to measure the vertical distribution of the soil temperature
and volumetric water content as function of time in the
lysimeter,
3. to measure ambient weather conditions to sufficiently
describe the boundary conditions for a coupled heat and
mass transfer model, and
4. to determine the physical characteristics of the soil
contained in the lysimeter.
Lysimeter Design. Installation and Calibration
Design and Construction
Weighing lysimeters were used in this research to meet the
objective of measuring hourly evaporation from a bare soil surface.
The weight of the lysimeter may be monitored by several different
methods. Mechanical measurement may be achieved by either supporting
the soil container with a lever and counterweight system or by
68
supporting the lysimeter directly by several load cells (Harrold, 1966;
Van Bavel and Meyers, 1961). The weight may also be detected by
supporting the soil container by a bladder filled with fluid and
monitoring the pressure of the fluid or the buoyancy of the lysimeter
container (Harrold, 1966; McMillan and Paul, 1961; King et al., 1956).
The hydraulic method has the disadvantages of being sensitive to
thermal expansion of the fluid and the high possibility of developing
leaks in the bladder and losing the fluid. The lever mechanism has the
advantage of requiring only one relatively low capacity load sensor
thus reducing instrumentation costs. The disadvantage of the counter
balance system is requirement of extensive underground construction to
contain the lever apparatus. This would require disturbing the
surrounding soil and border area and may be considered undesirable in
some cases. The cost of additional supporting structures and
excavation could be considerable as well. Three to four high capacity
load cells are required to support the lysimeter for direct weighing.
The higher cost of the load cells for the direct weigh method may
offset the cost of the increased excavation required for the counter
weight system. Directweigh lysimeters can be constructed such that a
high sensitivity can be achieved. Dugas et al. (1985) reported a
resolution of 0.02 mm of water for a direct weighing lysimeter with a
surface area of 3 m2. Most of the designs presented in the literature
have load cells which are difficult to access in the event of need for
maintenance or replacement.
A directweigh system was chosen because of its relatively simple
support system design and due to space limitations at the proposed
construction site at the Irrigation Research and Education Park (IREP)
located on the University of Florida campus in Gainesville, Florida
(Figure 31).
Aboukhaled et al. (1982) presented a relatively comprehensive
review of the existing literature regarding lysimeter construction and
design and noted the practical aspects of certain design
considerations. It was noted that the lysimeter annulus (gap area
between the soil container plus the area of the retaining wall) should
be minimized to reduce the effects of the thermal exchange between the
lysimeter soil mass and the surrounding air. Several researchers
(Aboukhaled et al., 1982; Black et al., 1968) made recommendations for
minimum lysimeter areas of 4 m2 for the primary purpose of reducing
the edge effect of the soiltosoil gap between the lysimeter interior
and the surrounding fetch. Maintaining a 100:1 ratio of fetch to crop
height has been recommended to minimize border effects. The site
selected for the lysimeter construction satisfied the fetch
requirements for very short crops or bare soil. Dugas et al. (1985)
noted that the problems caused by soiltosoil discontinuity should be
minimal if the lysimeter surface area were greater than 1 m2. Soil
depth should be sufficient so as not to impede root growth of crops to
be planted in the lysimeter.
Two soil containers were constructed using 4.8 mm (3/16 in.) steel
plate for the floor and side walls (Butts, 1985). All seams were
welded continuously to prevent water leaks and the entire container was
painted with an epoxybased paint to prevent corrosion. The lysimeter
measured 305 cm long, 224 cm wide and 130 cm deep. This provided a
soil surface area of 6.8 m2 and a soil volume of 8.9 m3. Main floor
supports consisted of four 20cm steel Ibeams (W8x15) and extended the
70
entire length past the edge of the soil container (Figure 32). To
alleviate the need for underground access to the load cells supporting
the lysimeter, two 2cm rods were inserted into both of the protruding
end of each floor joist and extended above the top of the lysimeter.
These rods were inserted through both webs of a 25cm steel Ibeam
(W10x26) extending parallel to the end of the lysimeter. The ends of
the rods were threaded such that nuts could be placed above and below
the upper Ibeam. The lower nut prevented the hanger beam (W10x26)
from sliding to the bottom prior to installation in the ground. The
top nuts prevented the rods from sliding out of the hanger beam after
installation. This assembly formed a cradle by which the lysimeter was
supported and provided access to four 4.5 Mg capacity load cells
located under each end of the hanger beams.
A pit with reinforced concrete block retaining walls and a 10cm
reinforced pouredinplace concrete floor was constructed for the
installation of the lysimeter at the IREP (Figure 32). The depth of
the pit was such that the bottom of the lysimeter tank was
approximately 15 cm above the pit floor and the tops of the retaining
wall and the lysimeter were flush. The pit floor was sloped toward the
center to drain water which might accumulate to a sump located
immediately outside of the North wall of the pit. Support columns were
incorporated into the retaining wall on which to place the load cells
supporting the lysimeter. The load cells could be installed or
removed by lifting the lysimeter with hydraulic automotive jacks from
ground level. Fill dirt was used to bring the surrounding soil up to
the elevation of the pit wall.
Prior to installing the lysimeters in the ground pits, two porous
71
ceramic stones (30 x 30 x 2.5 cm) were installed in the bottom of each
lysimeter and connected to vacuum tubing for the purpose of removing
excess water from the lysimeter. The lysimeters were then placed in
the pits and leveled using the top nuts of the supporting rods to
adjust the length of the rods and to maintain a fairly uniform load on
all of the rods.
The lysimeters were filled to a depth of 3 cm with a coarse sand
to provide free drainage of water to the ceramic stones. A Millhopper
fine sand excavated from a nearby site, was used to completely fill the
lysimeters. This soil was chosen because the profile was of fairly
uniform composition to a depth of 2 m and was typical of the local
sandy soils found in northern Florida and southern Georgia. After
loading the soil into the lysimeter, the soil was saturated with water
and allowed to stand overnight then removed by the vacuum system. This
provided settling to a density which might naturally occur in the field
after a long period of time.
Load cells manufactured by Hottinger Baldwin Measurements (HBM
Model USB10K) were then installed using a ball and socket connection
between the hanger beam and the load cell. A thin coat of white
lithium grease was applied to the surface of the ball and socket to
prevent corrosion of the contact surfaces. The output of the load
cells was a nominal 3 mV/V at full scale and could accept a maximum
supply voltage of 18 VDC. A regulated power supply provided 10 VDC
excitation voltage to each of the four load cells. Supply and output
voltages for each of the load cells was monitored using a single
channel digital voltmeter (Fluke Model 4520) and a multiple channel
multiplexer (Fluke Model 4506). Output was normalized using the ratio
72
of output voltage to input voltage (mV/V) to account for any variation
in supply voltage between load cells and time variation in the supply
voltage.
Lysimeter Calibration
Load cell calibration was performed by the manufacturer at the
factory for full load output of the cells. The manufacturer provided
load cell output (mV/V) at full scale load and no load. The load
carried by each load cell was expected to be near the full capacity of
each load cell with relatively small weight changes around that
reading. It was also noted that the calibration information provided
by the manufacturer was obtained under controlled atmospheric
conditions and using the factoryinstalled 3.3 m leads for measuring
the load cell output. The temperature sensitivity of the load cells
was very small compared to the changes in load cell output expected.
According the manufacturer's specifications, the effect of change in
load cell temperature upon the output was 0.08 percent of the load per
550C change in cell temperature. The installation site for the
lysimeters required the use of leads ranging in length from 37 to 53 m.
Therefore, an insitu calibration was necessary.
The soil surface of the lysimeter was covered with a polyethylene
sheet to prevent weight loss due to evaporation during the calibration
procedure. A single weight having a mass of approximately 12 kg was
added to the lysimeter and the normalized output of each load cell was
recorded. A second weight was added, thus increasing the total weight
added to the system. The load cell output was again recorded. This
was repeated until the cumulative weight added to the lysimeter reached
175 kg (equivalent to 26 mm of water distributed over the lysimeter
73
surface). The normalized output was then recorded as the weights were
removed.
To account for unequal distribution of weight due to variations in
soil density and slope of the lysimeter, the above procedure was
performed such that weights were placed directly over a single load
cell and repeated for each load cell and by placing the weight in the
center of the lysimeter as well.
Soon after the load cells were installed and calibration
completed, one load cell from each lysimeter was damaged by lightning
and removed for repair. Since during the calibration procedure, the
normalized output of the individual load cells was recorded,
regressions for the remaining load cells for each lysimeter could be
determined. The calibration was repeated after the load cells were
repaired and installed. Regressions for each of the lysimeters were
obtained for the calibration data with three and four load cells
(Table 31). It was determined in both cases that the slopes were not
significantly different at the 99 % confidence level for the two
lysimeters and thereby allowing the same regression for both lysimeters
to be used. Since, interest was in changes in weight and not in
absolute weight, the intercept would not be significant. The
calibration curves used in data analysis are shown in Figure 33.
Resolution of the lysimeters based upon the specifications of the
voltmeter, and standard error of the regression parameters was
calculated to be approximately 0.02 mm of water.
Temperature Measurement
The vertical distribution of temperature within the soil was
measured by monitoring the output of ANSI Type T thermocouples (copper
74
constantan) placed throughout the soil profile. Probes were
constructed using 2cm PVC pipe with small holes drilled radially into
the pipe at specified intervals (Table 32). Thermocouple wire
(22 AWG) was fed into the top end of the PVC probe and the end of the
thermocouple protruding through the radial holes. Epoxy cement was
used to secure the thermocouples in place. After the thermocouples
were installed in the probe, the lower end of the probe was cut at an
angle then plugged with wood and epoxy. Two identical probes were
installed symmetrically about the center of the lysimeter to provide
some degree of redundancy and to obtain some idea of the uniformity of
the temperature distribution at a given depth in the soil. Shielded
multipair thermocouple extension cable was used to carry the
thermocouple signal to the Fluke multiplexer and digital voltmeter.
Soil Water Content Measurement
The measurement of the distribution of soil moisture was needed to
determine initial conditions for the model and to verify the subsequent
simulation of water movement within the soil. Several methods were
considered for the measurement of soil water content within the
weighing lysimeters. The method chosen had to provide a relatively
accurate measure of the volumetric soil water content as a function of
depth with minimal disturbance of the surrounding soil. The chosen
technique must also provide data for a relatively small vertical soil
volume and should be commercially available.
Soil water content is very difficult to measure accurately
especially when trying to do so with minimal disturbance of the
surrounding soil. The simplest method for determining soil water
content is the gravimetric method. Core samples are obtained from
75
various depths, and placed in a water tight container, weighed and
placed in a drying oven. The oven temperature should be maintained at a
constant temperature (1050C) and the sample remain in the oven until a
constant weight is obtained. The sample is then removed from the
oven, allowed to cool in a desiccator, then weighed. This information
provides the water content on a weight basis. If the dry bulk density
of the soil is known, the volumetric water content can then be
calculated. The gravimetric method is the method by which all other
methods are calibrated (Schwab et al., 1966). However, the sampling
procedure disturbs the surrounding soil and the flow of moisture until
the soil sample is returned to the sample area. Even then the oven
dried soil is not at the same conditions as the undisturbed soil.
Therefore, the gravimetric method is not an acceptable method for
continuous or regular monitoring of the soil water profile.
Other methods for measuring the soil water content are classified
as indirect methods because they measure some property or
characteristic of the soil which varies predominantly with water
content. One of the simplest is the tensiometer. This instrument is
constructed of a hollow porous ceramic cup connected to a Plexiglas
tube. The tube is filled with water, sealed and connected to a
pressure sensor, usually a manometer, vacuum gauge or vacuum
transducer. After all of the air has been removed from the column of
water the tensiometer is placed in the soil such that the porous cup is
located at the desired depth. The water potential of the soil is
generally less than the potential of the water within the tube and is
drawn out of the tube until the water in the tube comes to an
equilibrium tension (or vacuum) with that in the soil. An electronic
76
vacuum transducer or manometer is used to measure the equilibrium
tension on the water column. The tensiometer has a range of
application from 0 to approximately 85 kPa (Long, 1982; Rice, 1969).
This would cover approximately ninety percent of the range of moisture
contents encountered in the field (Schwab et al., 1966). Using an
electronic pressure transducer would allow for the tension to be
monitored continuously. Several disadvantages in using tensiometers
exist. The soil water release curve (tension vs. soil water content)
must be known for each soil in which the tensiometer is installed.
However, this information is needed for the determination of the
concentration of water vapor within the soil for modeling purposes as
well. A second difficulty is that the tensiometer is essentially
detecting water content at a single point within the soil profile, thus
requiring the installation of several tensiometers within a finite
volume of soil. Tensiometers also require frequent maintenance. In
situ calibration of tensiometers is not required since the soil water
release curve is a unique property of the soil allowing the researcher
to obtain a soil sample from the experimental site and determine the
water release curve in the laboratory.
Electrical resistance is another property which varies with soil
water content. Two electrodes are imbedded within a porus gypsum block
which, when placed in the soil reaches an equilibrium moisture content
with the soil. The resistance between the two electrodes is measured.
However, as the water evaporates from the gypsum block, dissolved salts
remain in the gypsum. This deposition of salts within the porous
material will become redissolved when rewetted and cause an incorrect
reading. Their range of applicability extends mainly over dry soil
77
conditions (100 to 1500 kPa tension).
Nuclear emission methods provide yet another means by which the
soil moisture can be measured indirectly. In this method, either the
transmission of gamma rays or backscattering of neutrons is measured.
Fast neutrons are emitted from a neutron source and as they react with
the hydrogen nuclei present within the soil, they lose part of their
energy becoming slow neutrons. A detector is then used to count the
rate at which slow neutrons are backscattered in the soil surrounding
the probe. The neutron probe requires a single probe containing both
the radiation source and the detector. Access is gained by
installation of a tube in the soil thus allowing for repeated readings
at the same location with minimal disturbance of the soil after the
initial installation. However; the volume of soil penetrated by the
neutrons increases as the soil dries out, thus increasing the volume
over which the subsequent water content represents. Individual
readings are influenced by density and by the amount of endogenous
hydrogen. Once a calibration curve has been developed for the local
conditions, volumetric water content can be determined to within 0.5 to
1.0 percent moisture by volume (Schwab et al., 1966). A separate
calibration curve for measurements near the surface may also be
required due to fast neutrons escaping from the soil surface and not
being reflected back to the detector. Tollner and Moss (1985) noted
low Rsquare values of 0.6 and 0.4 for the calibration curves for a
neutron probe at a depths of 46 cm and 20 cm, respectively.
The use of gamma attenuation methods has gained some popularity
in recent years due to the ability of the instrument to measure the
volumetric water content or bulk density of relatively small volumes of
78
soil. The gamma probe requires the installation of two parallel tubes
spaced approximately 30 cm apart. A radioactive source which emits
moderately high energy gamma rays (Cs 37) is placed in one tube while a
detector is placed in the second tube at the same depth. The intensity
of gamma radiation reaching the detector is monitored by the detector.
Ferraz and Mansell (1979) discuss the gamma radiation theory used in
measuring soil water content in great detail. Primarily, the use of
the gamma attenuation apparatus in which the gamma rays are collimated
so that the intensity of the beam is focused and a high intensity
source are used. The use of uncollimated radiation in the field
requires the use of lower intensity sources due to the increased
possibility of radiation exposure to the operator. Ferraz and Mansell
(1979) state that the gamma attenuation method may be used insitu to
determine changes in volumetric water content quite accurately, but the
error associated with absolute moisture content can be quite high.
Ayers and Bowen (1985) used a gamma density probe designed for field
use to measure density of soil within a soil density box. Resolutions
of 0.016 g/cm3 for the wet bulk density were achieved in their
experiments.
The gamma attenuation method has several advantages over the
neutron method for measuring soil moisture. The volume of soil in
which the measurement is made remains constant regardless of water
content (Ferraz and Mansell, 1979) and consists of a band between the
source and detector probes approximately 5 cm deep (Ayers and Bowen,
1985). The disadvantage of the gamma method with respect to the
neutron method is the fact that a pair of tubes must be installed
instead of a single tube. However, the probes for the neutron and
gamma meters are typically the same diameter, thus allowing the use of
the neutron probe in the gamma access tube if necessary. The two
radiation methods will generally disturb less soil during installation
and will allow for recurring measurement at more frequent depth
intervals than with the tensiometers. Tensiometers also require
considerably more maintenance than do the radiation methods. The
tensiometers on the other hand can be monitored continuously by
automated data acquisition equipment while the radiation techniques
cannot. Extensive calibration is required for either of the radiation
measurement techniques. Based upon the desire to cause minimal
disruption of the soil profile and measure soil moisture at many points
within the soil, it was decided that a gamma attenuation technique
would be used to measure soil water content.
Gamma Probe Calibration
Extensive calibration curves were required to determine absolute
values of the soil water content using the gamma ray attenuation
method. While the soil was being loaded into the lysimeters, three
pairs of aluminum access tubes (5.1 cm O.D.) were installed in each
lysimeter (Figure 32). The tubes extended from the lysimeter floor to
23 cm above the soil surface. Parallel tube guides provided with the
Troxler 2376 Dual Probe Density Gauge were used to maintain parallelism
between the tubes as the lysimeter was filled.
The gamma probe primarily measures density of the test material;
therefore, it was necessary to measure the wet bulk density of the soil
corresponding to the depth at which the gamma readings were obtained.
A bulk density core sampler was used to obtain the soil samples for the
measurement of density and volumetric water content. The sampler
80
consisted of a series of removable brass rings of known volume within a
hollow cylinder. The cutting edge was tapered such that compaction of
the soil sample was minimized (Baver et al., 1972). The soil sample
was removed from the inner brass ring and placed in an aluminum soil
sample container and the lid closed to prevent moisture loss until the
samples from a single probe could be measured. Core samples were taken
at depth intervals of 10 cm from the soil surface to a depth of 71 cm.
Three probes were obtained from within the lysimeter to roughly
correspond to the locations of the paired access tubes. The soil
samples were weighed, dried, and reweighed. Wet density, dry density,
and volumetric water content were calculated for each sample. The soil
properties determined from each of the three probes were averaged for
each depth to account for the fact that the density between the access
tubes could not be measured. This procedure also provided information
regarding the areal uniformity of the soil with the lysimeter. The
dried soil was placed back in the holes in reverse order of their
removal so as to maintain the original soil profile. Care was taken
during subsequent sampling to avoid the soil core sites previously
sampled in the lysimeters.
The gamma readings were obtained in the following manner. The
meter was turned on for a minimum of twenty minutes prior to readings
being taken to allow for the circuits to stabilize. After warmup, the
detector probe and Cs137 gamma source were placed in the tubes of a
calibration stand. The calibration stand consisted of two parallel
aluminum tubes with materials of known density. The standard materials
were polyethylene (1.06 g/cm3), magnesium (1.75 g/cm3), magnesium
aluminum alloy (2.16 g/cm3), and aluminum (2.61 g/cm3). The source and
81
detector were placed in the stand level with the center of the
polyethylene standard and the amplifier gain of the meter was set so
that the peak count rate was achieved. The probes were then lowered to
correspond with the magnesium and a 4minute count was taken. In the
calibrate mode of the timer, the meter indicates the count per minute
for the 4minute interval. This was referred to as the standard count
and was used to normalize subsequent counts in the soil. The standard
count also provided a means by which to compare readings taken at
different times and accounted for variation in temperature of the gauge
and gain settings.
After the standard count was taken, the meter timer was set to
accumulate counts for one minute. The probes were placed in a pair of
tubes such that the source and detector were located at a depth
corresponding to the center of the core samples taken. A minimum of
two oneminute counts were recorded for each depth then the probe moved
to the next depth. After the last measurement at a depth of 66 cm was
obtained, the source and detector rods were removed from the access
tubes. Counts were taken at the same depths in all three set of access
tubes and the entire sampling process repeated in the second lysimeter.
Several attempts at this calibration procedure were made. Early
trials produced unacceptable ranges of scatter in the gamma probe
readings. At a later date, faulty electronic components were found and
expected to be the reason for the highly variable gamma probe readings.
After repair of the unit, the calibration procedure was repeated
yielding acceptable results.
According to Ferraz and Mansell (1979) the density of the test
material varies inversely with logarithm of the ratio of the intensity
82
of radiation observed in the test sample to the radiation intensity
measured in air. This relationship can be expressed as the count ratio
varying as a decreasing exponential with respect to the density of the
material. The count ratio is defined by the manufacturer as the ratio
of counts per minute observed for the test material to the counts per
minute observed in the magnesium standard (Troxler,1972). The observed
density of the soil obtained from the gamma measurement represents the
wet bulk density (Equation 31).
Pwet = A B ln(CR) (31)
where:
Pwet : wet bulk density of the soil [gcm3
A, B : constants of regression
CR : count ratio dimensionlesss]
count per min in test material
count per min in standard
If the dry bulk density is known, then the volumetric water content
can be calculated by
0 ( Pwet" dry) 32
Pw
where:
8 : volumetric water content [cm3cm3]
3
Pwet : wet bulk density of soil [gcm3]
Pdry : dry bulk density of soil [gcm3]
83
P : density ofyater [gcm3]
= 1.00 [g.cm ]
Combining equations (31) and (32) yields a relationship for
volumetric water content in terms of the count ratio and the dry bulk
density of the soil.
6 = A + Bln(CR) + CPdry (33)
It was conceivable that a calibration curve might be necessary for
each of the lysimeters and for various ranges of depth. Therefore,
linear regressions of the forms in Equations (31) and (33) were
determined for each lysimeter and for each of the depths. Statistical
testing revealed that a single regression could be used for all depths
and for both lysimeters for both the wet density and volumetric water
content. The final form of the calibration equation for determination
of volumetric water content was
0 = 1.52 0.331n(CR) 0.89Pdry (34)
The coefficient of variation (R2) for the regression shown in
equation (34) was 0.888. Measured volumetric water content when
plotted against the water content estimated by the calibration
(equation 34) should lie about a line with a slope of unity and an
intercept of zero (Figure 34). Examination of Figure (34) indicated
that relatively large errors could be associated with the exact
estimate of volumetric water content. Therefore, the gamma probe
84
should be used to determine changes in volumetric water content during
the experiments and gravimetric sampling should be used to determine
initial conditions for the experiment.
Soil Bulk Density and Porosity
The calibration procedure for the gamma probe provided sufficient
data to determine the vertical and horizontal distribution of the soil
dry bulk density in each of the lysimeters. As mentioned previously,
repeated core samples were obtained at several depths dispersed over
the area of each of the lysimeters. Figure (35) shows the average dry
bulk density as a function of depth in each of the lysimeters. The
error bars represent the standard deviation of the measured bulk
density. The bulk density for both lysimeters was relatively constant
at 1.38 g/cm3 in the top 25 cm and increased to approximately 1.6 g/cm3
between 25 and 40 cm. The largest deviation in measured bulk density
occurred at the lower depths.
The literature (Baver, 1972; DeVries, 1975) indicates that the
porosity for sandy soils ranges from 0.50 to 0.55 cm3/cm3. The soil
porosity was a soil parameter required as a function of depth. Due to
the variation in soil density, especially at the depths from 25 to
40 cm, it was decided that the porosity for the soil in the lysimeters
should be determined as a function of depth.
The porosity is defined as the volume of pore space in the soil
per unit of bulk volume. An air pycnometer (Baver, 1972) was used to
measure the porosity of the soil. Four levels of soil bulk density
ranging from 1.2 to 1.6 g/cm3 were used for the tests. Three soil
samples at each level of bulk density were prepared from a soil sample
obtained from the weighing lysimeters. The volumetric water content
85
of each soil sample was 0.05 cm3/cm3. The mass of soil placed in the
50 cm3 sample cup of the air pycnometer was that required to achieve
the prescribed levels of dry bulk density. The test was performed for
the sample in the cup and the void volume recorded. The volume
occupied by the water in the sample was added to the measured void
volume to yield the total pore space in the soil. A total of three
measurements were made for each sample. The test was also repeated
for a single soil sample with the maximum density achievable
(1.64 g/cm3). The porosity was calculated for each measurement by
dividing the 50 cm3 sample volume into the void volume and were
averaged for each bulk density level (Table 33). The repeatability of
the experiment was indicated by the standard deviation of the porosity
for each of the samples, while the accuracy of the measurement was
indicated by the standard deviation for the average of all measurements
for a given density level. The porosity was found to vary linearly
with bulk density (Figure 36). The error bars shown in Figure 36
indicate the standard deviation of the soil porosity at each bulk
density level. The porosity measurements were then used in defining
the soil profile for the validation and calibration simulations.
Experimental Procedure
Experiments were designed to monitor the energy and mass transport
processes from a bare soil surface for two weighing lysimeters on a
continuous basis. Automatic data collection was performed by a Digital
Electronics Corporation (DEC) PDP11/23 minicomputer equipped with an
IEEE488 interface board and utilizing the RT11 operating system. The
RT11 operating system provided several system subroutines which were
primarily for scheduling the data collection as well as other time
86
manipulations. The main program and associated subroutines were
written in FORTRAN66. Subroutines for initializing and retrieving
data through the IEEE interface were written in either assembly or
FORTRAN by Dr. J. W. Mishoe for previous research. All signals from
the sensors used were analog signals and were monitored on 50 channels
of a Fluke 4506 multiplexer and subsequently read via the IEEE
interface from a Fluke 4520 digital voltmeter. Data monitored by the
computer for each lysimeter were net radiation, soil temperature
distribution, and supply and output signals from each of the load
cells.
Net radiation was measured using net radiometers (WEATHERtronics,
Model 3035) mounted 85 cm above the soil surface of the lysimeter. The
net radiometers utilized blackened thermopiles as the sensing element
to detect the difference between incoming and outgoing radiation. A
positive millivolt signal was produced when the incident radiation was
greater than that being reradiated, while a negative signal indicated
that more energy was being radiated from the surface than impinging
upon the surface. Factory calibration curves were used to convert the
millivolt signals to net heat flux (W/m2). The millivolt signal was
read at two minute intervals and total millivolts and number of times
read were recorded on floppy diskette at ten minute intervals.
Thermocouples, supply voltage, and load cell output were recorded at
ten minute intervals as well.
Data files were closed on an hourly basis and new ones opened for
each hour. This minimized the risk of data loss due to power outages
or other computer malfunctions. The most data that would be lost that
had already been recorded would be that for one hour in the event that
the power was lost just prior to closing the data file. When power was
restored, the PDP11/23 would automatically restart and begin the data
collection routine. Separate data files were maintained for each of
the two lysimeters.
Hourly values of the ambient dry bulb temperature and relative
humidity were measured within a standard weather shelter from the
adjacent weather station monitored and maintained by the Agricultural
Engineering Department. Sensors for measuring air temperature were
Type T thermocouples and instantaneous values were recorded hourly. A
Campbell Scientific CR207 sensor was used to monitor relative
humidity and consists of a wafer whose electrical properties vary with
relative humidity. Manufacturer's literature states that the sensor
provides reliable information when used in noncondensing conditions
with a relative humidity between approximately 20 and 90 percent.
However, the wafer material tends to absorb moisture over time
decreasing the reliability of data obtained. To account for the lag
time of the sensor, data was averaged over an hour and recorded.
Relative humidity data was generated as well by assuming that the
minimum daily temperature was the dewpoint temperature for the day.
Linear interpolation between consecutive daily minimums to provide a
continuous estimate of the dewpoint temperature throughout the day.
Relative humidity could then be calculated based upon the estimated
dewpoint temperatures. Wind speed and direction were measured at a
height of 2 m within the same weather station. Meteorological data
were recorded hourly and uploaded daily to the VAX mainframe managed by
IFAS. Access to the hourly data was achieved through the AWARDS
(Agricultural Weather Acquisition Retrieval Delivery System).
Meteorological data, including relative humidity determined from the
estimated dewpoint temperature, are shown in APPENDIX B.
Tests were begun by an irrigation or rainfall event to ensure that
the upper layers of the soil were sufficiently wet so as to provide a
minimum of three to four days of evaporation data. Most frequently,
water was added to the soil the evening prior to beginning the test to
insure that the soil surface was wet and to allow some downward
distribution of water prior to the initiation of each test. Core
samples were obtained to determine the initial distribution of water
within the lysimeter. Samples were taken at vertical intervals of
5.1 cm for the first 30.5 cm then every 10.2 cm until a depth of
91.4 cm was reached. Gamma probe measurements were taken in each of
the three pairs of access tubes corresponding to the center of each of
the core samples. Actual depths of gamma probe measurements and core
samples are shown in Table 34. Subsequent gamma probe measurements
were made approximately every other day. More frequent readings were
not obtained in most cases due to the time required for each set of
readings (2 hrs.). This minimized disruption of the other data being
collected.
Obtaining the gamma probe readings required placing a scaffold
across the lysimeter to avoid disturbing the weight of the lysimeter
and compacting the soil. This, in turn, could have disrupted the
evaporation process due to shading of the soil surface as well as some
of the net radiation measurements. These periods when the integrity of
some of the data may have been questionable was marked in the hourly
data file by a flag entered via the keyboard. The flag was turned on
when gamma readings were being made on each lysimeter and turned off
again after the task was completed.
Data were collected continuously over a period extending from
November, 1986 to August, 1987. Data collection was interrupted on
several occasions due to equipment failures, power outages, equipment
maintenance, and instrument calibration. Most tests were three to four
days in length with a few extending to six or seven days.
Data Analysis
Data files containing the data recorded at tenminute intervals
were concatenated into daily files. The daily files were then combined
to correspond to specific dates for the individual tests. Data
contained in the files consisted of dayofyear, timeofday, the
normalized output for each load cell (mV/V), the 10minute total net
radiometer output, the number of net radiometer readings during the
10minute interval, the status flag, and twenty temperatures
corresponding to the depths indicated in Table 32.
The normalized output of the load cells was totaled and the
average total millivolt per volt for the hour was calculated. The
difference between successive hourly totals was used in the load cell
calibration equations to determine the hourly and cumulative
evaporation of water from the lysimeter. The output signal for the net
radiometers was totalled for the hour then an average rate of net heat
flux was determined using the manufacturer's calibration curves. The
two temperatures for each depth were averaged over the hour as well.
If the data flag had been turned on at any time during the hour, it was
assumed to have been on the entire hour. Reduced data files contained
the hourly values of julian date, time, cumulative evaporation (mm),
hourly evaporation (mm/hr), net radiation (W/m2), data flag, and
