Title Page
 Table of Contents
 List of Figures
 Energy and water transfer in a...
 Equipment and procedures for measuring...
 Measuring thermal diffusivity of...
 Model analysis
 Summary and conclusions
 Biographical sketch

Title: Modeling the evaporation and temperature distribution of a soil profile
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00089826/00001
 Material Information
Title: Modeling the evaporation and temperature distribution of a soil profile
Series Title: Modeling the evaporation and temperature distribution of a soil profile
Physical Description: Book
Language: English
Creator: Butts, Christopher Lloyd
Publisher: Christopher Lloyd Butts
Publication Date: 1988
 Record Information
Bibliographic ID: UF00089826
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001125835
oclc - 20111916

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
    List of Figures
        Page v
        Page vi
        Page vii
        Page viii
        Page ix
        Page x
        Page xi
        Page 1
        Page 2
        Page 3
        Page 4
    Energy and water transfer in a sandy soil: model development
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
    Equipment and procedures for measuring the mass and energy transfer processes in the soil
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
    Measuring thermal diffusivity of soils
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
    Model analysis
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
        Page 197
        Page 198
        Page 199
        Page 200
        Page 201
        Page 202
        Page 203
        Page 204
        Page 205
        Page 206
        Page 207
    Summary and conclusions
        Page 208
        Page 209
        Page 210
        Page 211
        Page 212
        Page 213
        Page 214
        Page 215
        Page 216
        Page 217
        Page 218
        Page 219
        Page 220
        Page 221
        Page 222
        Page 223
        Page 224
        Page 225
        Page 226
        Page 227
        Page 228
        Page 229
        Page 230
        Page 231
        Page 232
        Page 233
        Page 234
        Page 235
        Page 236
        Page 237
        Page 238
        Page 239
        Page 240
    Biographical sketch
        Page 241
        Page 242
        Page 243
Full Text







112- S


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 well-rounded education of

my fellow graduate students, particularly, Bob Romero, Ken Stone, Matt

Smith, Ashim Datta and Kumar Nagarajan.



ACKNOWLEDGEMENTS . . . . . . . . . . . . ii


ABSTRACT . . . ...... . ....... . . . . . .x


I. INTRODUCTION . . . . . . . . . . . 1

Problem Statement . . . . . . . . 1
Research goals and objectives . . . . . . 3


Introduction . . . . . .
Literature Review . . . .
Model Objectives . . . . .
Model development . . . .
Determination of Model Parameters
Numerical solution . . . .



Introduction . . . . .
Objectives . . . . . .
Lysimeter Design, Installation a
and Construction . . .
Experimental Procedure . . .
Data Analysis . . . .
Experimental Results . . .
Summary . . . . . .

. . . . . . .o

. . . . . . .

d Calibration Design

. . . . . . .
. . . . ... .
. . . . . . .
. . . .
. . . . .


Introduction . . . . . . .. . . .
Objectives . . . . . . .. ....
Literature Review . . . . . . . . .
Procedure . . . . . . . . . . .
Results and Discussion . . . . . . ...
Conclusions . . . . . . . . . . .







. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

V. MODEL ANALYSIS . . . . . . . . .. . .


Introduction . . . . . . . . .. . 158
Validation . . . . . ...158
Sensitivity Analysis . ..... . . . .. 170
Summary . . . . . . . . . . 176

VI. SUMMARY AND CONCLUSIONS . . . . . . . . . 208

BIBLIOGRAPHY . . . . . . . . .. . . . . 212

TRANSFER MODEL . . . . . . . . . . . 219

STUDIES . . . . . . .... . . . . . 225

BIOGRAPHICAL SKETCH . . . . . . . . . . . 241

Variable Definition
B empirical function for determination of Cer and Chr
Cdr surface drag coefficient dimensionlesss]
Cer Dalton number, dimensionless mass transfer
Chr Stanton number, dimensionless heat transfer

Cpa specific heat of moist air [J'kg-1"K-1]
Cps specific heat of solid constituent of soil mixture

Cpw specific heat of water [J-kg-1'K-1]
Cpv specific heat of water vapor [J-kg-1-K-1]
Cs volumetric heat capacity of soil mixture (solid,
liquid and air phases) [J'kg-l1(-1]
Da diffusivity of water vapor in air [m2"s-1]
DL unsaturated hydraulic diffusivity of soil [m2.s-1]
Dt thermal diffusivity [m2.s-1]
Dv diffusivity of water vapor in soil [m2.s-1]
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
[W-m- ]

g gravitational acceleration [m's-2]

gj shape factor in the j-th principal axis for thermal
conductivity calculations dimensionlesss]

H vertical component of heat flux in the atmosphere
[J-m ]

hfg latent heat of vaporization of water [J3kg-1]
hh convection heat transfer coefficient for soil
surface [W -m'2K-i]

hm convection vapor transfer coefficient for soil
surface [ms- ]

K hydraulic conductivity [m0'm-2.s-1]

Kh, Km, Kw turbulent transfer coefficients for peat, momentum
and water vapor, respectively [m s ]

kij ratio of temperature gradient in i-th soil
constituent to the temperature gradient in the
continuous constituent (water or air) in the
direction of the j-th principal axis dimensionlesss]

LE latent energy transfer in the atmosphere [Jm-2]
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 [kg-m- '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 [W-m- ]

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'm-2]

qsj,n rate of sensible heat change of node j at time step
n [W'm-3]

qvj+1,n rate of heat carried into node j from node j+1 by
water vapor movement at time step n [W'm-2]
R universal gas constant [kg-m2-s-2.mol-1(-o-1]

Rw gas constant for water vapor, determined by dividing
the universal gas constant, R, by the molecular
weight of water [m2.s-2OK-1]
Re Reynolds number dimensionlesss]
Ri Richardson number dimensionlesss]
Rn net r diation flux incident upon soil surface
[W-m- ]
S soil porosity [m3m-3]

s slope of the saturated vapor pressure line

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 [ms-l]
U* shear velocity [m-s-1]

V volume [m3]

WS Wind speed [mss-1]

xi volume fraction of the i-th soil component [m3m-3]
z depth below soil surface [m]
z0 surface roughness length [m]
a soil tortuosity dimensionlesss]

I- psychrometric constant [kPaCK-1]
7r diabetic influence function dimensionlesss]

soil water potential [m]
X thermal conductivity of soil composite at node j
v kinematic viscosity [m2.s-1]

Pa density of moist air [kg'm-3]

Ps dry bulk density of soil [kg'm-3]
Pv water vapor concentration [kgH20'ms3il]

Pvj,n water va or concentration at node j and time step n
Pva water va or concentration of ambient air

Pvs water vapor concentration in the air at saturation
Pw density of water [kg-m-3]
9 dimensionless water content


Variable Definition (continued)
0 volumetric water content of soil [m3.m-3]

Ojn volumetric water content of soil at node j and at
time step n [mn3m-3]

Or residual volumetric water content [m3.m-3]
es volumetric water content at saturation [m3m-3]
7 shear stress [N-m-2]

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



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.


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



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, no-till and row


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


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

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.



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



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


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 Blaney-Criddle

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 Blaney-Criddle 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


cpaPa (e(Ts) e(Td) )
E = Hf- R (2-1)


E = evaporation rate [m-s-1]

Cpa = specific heat of air [J-kg-l'o-1]

Pa = density of dry air [kg'm-3]
hfg = latent heat of vaporization [J'kg-1]

= psychrometric constant [kPa*'K-1]

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


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 well-watered 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 soil-atmosphere 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


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,


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 (2-2)

S = p aKm aU (2-3)

H = cpaaKh aT (2-4)


Qv = vertical flux density of water vapor [kg-m-2.s-1

Pa = density of air [kg'm-3]
7 = shear stress rate [N'm-2]

H = vertical flux of heat [J.m-2]

z = vertical distance [m]

q = specific humidity [kgapor/kgdry air]
U = velocity [m's-1]

Cpa = specific heat of dry air [J*kg-1K-l]

T = air temperature [OK]

Kh, Km, Kw = turbulent transfer coefficients for Peat, momentum
and water vapor, respectively [m' s- ]

Combining the equations 2-2 and 2-3 and rearranging the following

expression for evaporation is obtained.

Q V Kwm ~ (2-5)

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

Monin-Obukov length (L) are used as a measure of atmospheric

instability and can be determined by equations 2-6 and 2-7,


-T- (-- +r)
Ri = ---aU- (2-6)

L (2-7)
( kgH

= acceleration due to gravity [m-s-2]

r = adiabatic lapse rate [(K(m-1]
S 9.86 x 10-3

k = von Karman's constant dimensionlesss]
= 0.428

U* = shear velocity [mns-1]

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 60-minute

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

7 = PaKm,- awu (2-8)

H = cpaPaKh + cpaPaw'' (2-9)

LE = -hfg(paKw paw'q ) (2-10)
T7 : time averaged turbulent fluctuation of the air
temperature [K(]
u7 : time averaged turbulent fluctuations pf the horizontal
component of the wind velocity [m-s"]
w : wind velocity in the vertical direction [m's-1]
-w : time averaged turbulent fluctuations of the vertical
component of the wind velocity [m-s-11

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 (2-11)

LE = hfgPaw (2-12)

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.


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 2-4 and 2-2), referred to as the Bowen ratio (BR), is given


H cpaPa Kh Jt
BR = (2-13)

The energy balance at the soil surface is

Rn S G LE = 0 (2-14)


Rn = net radiation incident upon soil surface [W*m-2]

S = soil heat flux [W-m-2]

G = sensible energy flux into atmosphere [W-m-2]

LE = flux of latent energy from soil [W-m-2]

An expression for the rate of evaporation of water from the soil can be

determined by substitution of equation (2-13) into (2-14) 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


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) (2-15)


Ep potential evaporation [m3.m-2]

ea = vapor pressure of atmosphere [kPa]

eas = saturated vapor pressure at the air temperature [kPa]
h = surface vapor transfer coefficient [m's-1]

G = sensible heat flux in the atmosphere [W-m-2]

Rn = net radiation upon the soil surface [W-m-2]
s = slope of saturated vapor pressure line [kPa.'K-1]

7 = psychrometric constant [kPaOK-1]
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


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


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 (2-14). 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.

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

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- (2-16)

aT a aT
Cs at = ( )+ Q(O,T) (2-17)

0 = volumetric water content [m3.m-3]
t = time [s]
z = depth below soil surface [m]
T = soil temperature [OK]
Cs = volumetriq heat capacity of soil mixture
S= thermal conductivity of soil mixture [W-m-1OK-I]
K = unsaturated hydraulic conductivity [m's-1]

Do = hydraulic diffusivity [m2.s-1]

Sao0 D a VS9Pw
S- + PRT a80
Da = diffusivity of water vapor in air [m2.s-1]
P = atmospheric pressure [kPa]
v = ratio of atmospheric pressure (P) to partial
pressure of air
S = soil porosity [m3m-3]

g = acceleration due to gravity [m-s-1]
= soil matric potential [m]
Rw = ideal gas constant for water vapor

Pw = density of water vapor [kg'm-3]


DT = diffusivity of water (liquid and vapor) due to
temperature gradients [m*'s-10'1]

+ Da vSahon
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'm-3"K-1]

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


Jackson (1964) studied the non-isothermal movement of water using

equations (2-16) and (2-17) 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


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


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

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


development of the model for coupled heat and mass transfer in the


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) (2-17)


Pv = mass concentration of water vapor in the soil air
space [kg-m-3]

Dv = Diffusivity of water vapor in the soil air space
= Da

Da = Diffusivity of water vapor in the air [m2.s-1]
a = soil tortuosity dimensionlesss]

Ev(z,t) = rate of phase change of water [kg'm-3-S-1]

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


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


at-W) = aQ EL(z,t) (2-18)


0 = volumetric water content [m3m-3]

Pw = density of water [kg-m-3]

QL = diffusion of liquid water [kg-m-2 s-]

EL = rate of change of liquid water to water vapor
[kg-msils- ]
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- (2-19)

qL = volumetric flow rate of liquid water [m3.m-2s-1]
k = unsaturated hydraulic conductivity [m-s-1]

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 (2-19) as follows

L -k a (2-20)
The hydraulic diffusivity of the soil is defined as

DL M K -a (2-21)
By substituting the definition for the hydraulic diffusivity (2-21) and
the volumetric flow rate in terms of the soil water content (2-20) the
equation for the conservation of mass for liquid water (2-18) becomes

aB a a( E(z,t) (2-22)
aT = az (DL -) Pw

The conservation equations for the liquid (2-22) and vapor (2-17)
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) (2-23)
Rearranging equation (2-23) and solving for Ev(z,t) yields

Ev(z,t) = Ez,t) (2-24)

Substitution of equation (2-24) into the vapor continuity equation
(2-17) results in the following equations describing the conservation
of water vapor within the soil.

SPv a aPv Ev(z,t)
= (D t ) (S- ) (2-25)

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 (2-26)

X = thermal conductivity of soil mixture [Wnm-I.O-I]

T = soil temperature [OK]

hfg = latent heat of vaporization [J-kg-1]
Cs = volumetric heat capacity of soil mixture [J-m-1K-1]
= (1-S)Pscps + Opwcpw + (S-O)Pacpa
S = Soil porosity [m3m-3]

ps = dry bulk density of soil [kg'm-3]

Pa = density of moist air [kg-m-3]
cps = specific heat of solid portion of soil [J-kg-loK-1]

Cpw = specific heat of water [J-kg-l.'K-1]

Cpa = specific heat of moist air [J-kg-1oK-1]
Equations (2-22), (2-25), and (2-26) 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

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


vapor density can be determined from the saturated vapor pressure using

the ideal gas law.

pV = nRT



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

R, T



Pvs = concentration of water vapor t saturated vapor pressure
per unit volume of air [kg-m" ]

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.s-2K-1]

The above equation (2-28) is for the case in which the chemical

potential of the water is zero. However, due to matric and osmotic


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(- -) (2-29)


ev(0) = saturated vapor pressure of water with chemical
potential, 0 [Pa]

= soil water potential [m]

g = acceleration due to gravity [mns-2]
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(- -) (2-30)

A set of simultaneous equations (2-22, 2-25, 2-26, 2-30) 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 (2-30) and provides a fourth state

equation to be used in the model.

A well-posed 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 one-dimensional and

therefore required two boundary conditions. The first boundary is the

obvious soil-atmosphere 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) (2-31a)

-DL --= P (2-31b)

-Dv = hm(pvO- Pva ) (2-31c)


cpw = specific heat of water [J-kg-1.''1]
Cpv = specific heat of water vapor [J-kg-1.K-1]

DL = hydraulic diffusivity [m2s-1]
Dv = diffusivity of water vapor in soil [m2s-1]
hh = boundary layer heat transfer coefficient [W-m-2.K-1]
hm = boundary layer mass transfer coefficient [m-s-l]
P = precipitation or irrigation rate [m3.m-2s-1]

Rn = net radiation incident upon soil surface [WNm-2]
T = soil temperature [oK]

Tavg = average temperature [OK]
X = thermal conductivity of soil [W.m-l1.- 1]

Pva = mass oJ water vapor per unit volume of dry ambient air

PvO = concentration of water vapor at the soil surface, z=0
0 = volume of water per unit total soil volume [m3.m-3]

Ideally, for a semi-infinite 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 semi-infinite 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 (2-32a)
az Zzo
o I = 0 (2-32b)
a zz0Z

z = 0 (2-32c)

The system defining the movement of water vapor, liquid water and

energy is defined by a system of partial differential equations (2-22,

2-25, 2-26) with boundary conditions (2-31a-c, 2-32a-c). 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

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 10-5 (2-33)
va p


Dva = diffusivity of water vapor in air [m2.s-1]
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 10-5 (2-34)
~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 (2-34).

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 (2-35)


C = volumetric heat capacity [J-m-3.K-1]

x = volume fraction [m3.m-3]

q = quartz

m = mineral

o = organic


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 2-1). 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 (2-36)
Sqxq + kmXm + koXo + kX + kaxa


S = thermal conductivity of the soil [W-m-1K-1]

\ = thermal conductivity of soil component [W'm-1.K-I1]
xi = volume fraction of soil component [m3m-3]

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 i-th 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) (2-37)

kij = the ratio of the temperature gradient in the i-th
component to the temperature gradient in the continuous
component in the direction of the j-th 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
The ratio of the temperature gradients in the i-th component can be
determined by the following:

kj = 1 + 1 (2-38)
+ ( 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 2-1). 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 (2-37) 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


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 (2-35) for

the volumetric heat capacity and equations (2-36), (2-37) and (2-38) 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 2-1).

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 2-2), while the specific water

capacity is the slope of the soil water retention curve. The

hydraulic conductivity (Figure 2-3) and the specific water capacity

(Figure 2-2) 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

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 (2-39)
Ts r 1 + (a)n


9 = dimensionless water content

0 = volumetric water content [m3m-3]

r = residual volumetric water content [m3m-3]

es = volumetric water content at saturation [m3m-3]
= 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) [(1-e1/m)-m+ (1_01/m)m_2] (2-40)

The Van Genuschten method yields a continuous function for the
hydraulic diffusivity which is highly desirable for numerical


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 potential-water 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 A-i and A-2 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 (2-41)


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


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) (2-42)


WS = wind speed [ m-s-1]

U* = shear velocity

= TI

z0 = surface roughness height [ m ]
k = von Karman constant dimensionlesss]

T = shear stress at the surface [N-m-21
p = density of air [kg-m-3]
According to Sutton (1953), for practical purposes, the shear velocity

can be estimated as

U = (2-43)

The bulk mass transfer within the surface sublayer is defined by

Qm = CeyWSr (Pvs Pvr) (2-44)


Qm = vertical mass flux [kg'm-2.s-1]

WSr = wind speed at reference height, zr [m's-1]
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'm-3]
The surface mass transfer coefficient (hm) used in this model is

related to Cer by

hm = WS, Cer (2-45)

Brutsaert (1982) presented functions for the Dalton number in terms of

the drag coefficient and the dimensionless Schmidt number as

Cer = (2-46)
(B + Cdr)


B = function of dimensionless Schmidt number

Cer = Dalton number

Cdr = surface drag coefficient dimensionlesss]




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 (2-47). The surface is rough if the Reynolds number is

greater than 2.0 and Equation (2-48) is used to determine the value

of B.

B = 13.6 Sc 13.5 (2-47)

/4 2/
B = 7.3 Re 4 Sc 5.0 (2-48)


Sc = Schmidt number

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 (2-47) and (2-48)

above. The surface heat transfer coefficient can then be determined by

hh = Chr Pacpa WSr



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

h = k2WS ( 7 + In (z + D) (2-50)
m Zo


hm = transfer coefficient of mass from surface to height, z,
in the air [m-s-1]

k = von Karman constant dimensionlesss]

WS = wind speed at height, z [m-s-1]

D = height displacement [m]
= d + z0

d = height above soil surface where wind velocity is zero
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


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 2-4). 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 (2-50) 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 2-50) and heat

(Equation 2-51) 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 2-2

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


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) (2-52)
ax xj+1 xj

af f(xj) f(xj-1)
backward: = Vf = (2-53)
Xj xj-1
central: = 6f (2-54)

1 f(xj+1) f(xj) f(xj) f(xj-1)
T xj+1 xj xj xj-1

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 2-5).
The partial differential equation defining the conservation of
energy within the soil profile (equation 2-26) 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


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.

dx-dydwj)qsj = dx-dy (qcj+l,n + qcj-l,n + qLj+l,n

+ qLj-l,n + qvj+l,n + qvj-l,n)



The subscripts

s :
c :
L :
j :

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


dx-dy :
dx-dyddwj :
dwj :

The change in

to n+1 is

cross-section 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


The heat flux conducted from nodes j+1 and j-1 to node j at a given

time step, n, is described by Fourier's law of heat conduction and is

expressed as follows

qcj+l,n= (Tj+i,n- Tj,n) (2-57)

qcj-l,n= (Tj-l,n- Tj,n) (2-58)

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 (2-59)

(Pvj+,n Pvj,n) (Tj+1,n+Tj,n)
qvj+1,n = cpvjDvj dz 2 (2-60)

and similarly for the heat flux transferred by mass diffusion from node
j-1 to node j as

(0j-l,n 1j,n) (Tj-l,n+Tj,n)
qLj-l,n = PwjCpwjDLj dzj-1 2 (2-61)

S(Pvj-l,n Pvj,n) (Tj-l,n+Tj,n)
qvj-l,n = cpvjDvj dzj- 2 (2-62)

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



qej,n =

Substitution of the numerical expressions for the individual terms of
the energy balance (equation 2-55) and dividing by the cross-sectional
area normal to the z axis (dx-dy) yields

Tj,n+l-Tj,n Tj+l,n-Tj,n Tj-1,n Tj,n
dwjCs dt= J dz- + dzj_

(Oj+1,n-0j,n) (Tj+l,n+Tj,n)
+ PwjCpwjDLj dzj 2

(0j-i,n 0j,n) (Tj-l,n+Tj,n)
+ PwjcpwjDLj dzj- 2

(Pvj+1,n Pvj,n) (Tj+1,n+Tj,n)
+ cpvjDvj dzj 2

(Pvj-l,n Pvj,n) (Tj-l,n+Tj,n)
+ cpvjDvj dzj-. 2

dwj (hfg,jEj,n) (2-64)

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 j-1 to node j less the amount of water changed to a
vapor phase. The numerical expression for the water continuity per
unit area becomes


(0j,n+l- Oj,n) 0 +1,n- 0j,n (Oj-l,n- Oj,n)
dwJ dt DLj dzj + dzj-l

dwj (2-65)

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
(Pvj-l,n Pvj,n) dwj Ej,n
+ -Dv dzj-1 + (Sj j,n)

Equations (2-64), (2-65), and (2-66) 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. 2-1), 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

water: --

vapor: -T-

(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



: net solar radiation incident
at time step n [W-m'2]

: precipitation at time step n
: boundary layer mass transfer
step n [ms' ]

upon the soil surface

coefficient at time

: boundary layer heat transfer coefficient at time
step n [W-m- 2.- ]
: air temperature at time step n [JK]
: ambient water vapor concentration at time step n

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 zero-flux condition

can be represented. These will be developed for the water for the

purpose of illustration.







If the water content of the last node (j=nc) were equal to that of
the preceding node (j=nc-1), 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(nc-l,n Onc,n) (2-70)

If the distance between nodes nc and nc-1 is the same as the distance
between the nodes nc and the imaginary node nc+1, then it follows that

Onc+l,n = Onc-1,n (2-71)

Substitution of 2-71 into 2-65 yields as the boundary condition for the

node, j=nc, the following

(Onc,n+1- Onc,n) DLnc
dwnc dt 2 d-- (Onc-1,n-nc,n)
dwnc wnc

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 (Pvnc-1-Pvnc)
dwnc Enc,n

(S Oj,n)

(Tnc,n+1-Tnc,n) (Tnc-l,n- Tnc,n)
dwnc Cs,nc dt = 2 fnc dznc-1

(nc-1,n Onc,n) (Tnc-l,n+ Tnc,n)
+ 2 PwncCpwncDLnc dznc-1

(Pvnc-l,n'Pvnc,n) (Tnc-l,n+ Tnc,n)
+ 2 CpvncOvnc dznc-- 2

+ dwnc Enc,n (2-74)

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 (2-64) through (2-69)and (2-72) through
(2-74) along with equation (2-30) 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

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 round-off

or truncation errors and would not propagate or grow with time and

would yield an unconditionally stable solution (Conte and de Boor,


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 (j-1) 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 j-1 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

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


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 2-40) and thermal properties (equations 2-35 and

2-36) of the soil and the surface transfer coefficients (equations 2-50

and 2-51) were calculated at the beginning of each time step. A

general description of the solution algorithm is shown in Figure (2-6).

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 (2-30). 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 2-3). The program was then

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 2-7 and 2-8, 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 2-1.


Formulae for the geometric shape factors used in
calculating thermal conductivity based upon ratios
of the unit vectors of the principal axes.


Ellipsoid of revolution
x=y=nz n = 0.1
n = 0.5
n = 1.0
n = 5.0
n = 10.

Elongated cylinder with
elliptical x-section

x = ny;

0.14 0.14



(n + 1)-1



n/((n +1)





Typical sand grain

Table 2-2. 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) 2-22

Conservation of mass (vapor) 2-25

Conservation of energy 2-26

Vapor-Liquid equilibrium 2-30

Surface Boundary Conditions

Conservation of mass (liquid) 2-31b

Conservation of mass (vapor) 2-31c

Conservation of energy 2-31a

Lower Boundary Conditions

Conservation of mass (liquid) 2-32b

Conservation of mass (vapor) 2-31c

Conservation of energy 2-32a

Soil properties

Diffusivity of water vapor in soil (Dv) 2-34

Volumetric heat capacity (Cs) 2-35

Thermal Conductivity (X) 2-36

Hydraulic Diffusivity (DL) 2-40

Surface Mass Transfer Coefficient (hm) 2-50

Surface Heat Transfer Coefficient (hh) 2-51

Table 2-3.

Mesh spacing used in grid generation for the coupled
heat and mass transfer model for sandy soils.

Cell Width



Distance to
next node












o 3


0.0 0.1 0.2 0.3 0.4 0.5 0.6

Volumetric Water Content (m3.m-3)

Figure 2-1. Thermal conductivity of a typical sandy soil as a function
of volumetric water content (DeVries, 1975).


104 0.010

1 \ _-- Capacity -/ 0.008 E
E 103 -
/ ooo8


0.004 )
o /
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 2-2. Soil water potential (tension) and specific water capacity
as a function of volumetric water content for a sandy soil.



10 -- Diffusivity 10
~1 Conductivity o105
E 104
1 0-- 103
10-2 ---
0 102

o 10-4..
o 10-
10-6-- /*10-2 1
1 10 ..10-3

10-8 I 10o-4
0.00 0.10 0.20 0.30 0.40
Volumetric Water Content ( cm3/cm3)

Figure 2-3. 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


Richardson Number ( Ri)

Figure 2-4. Relationship of the diabetic influence function and the
Richardson number.





I i ,


-1.0 -

P Rn Qv=f(To,RH,hm,WS)



AZ j- j-l




Figure 2-5. Schematic of discretized domain for evaporation and soil
temperature model.

Figure 2-6. Flowchart describing the solution algorithm to solve the
system of finite difference equations for a coupled heat
and mass transfer model.



o" 1



60 120
Time Step ( s )



Figure 2-7. Sum of squared residuals for the soil profile energy
balance after a 72 hour simulation using various time


_1_1111 _



60 120
Time Step (




Figure 2-8. Sum of squared residuals for the soil profile mass balance
after a 72 hour simulation using various time steps.




if fll





The process of evaporation of water from the soil is a complicated

process in which the dynamic processes of heat and mass transfer are

inter-related. 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 non-weighing. As their names imply, weighing lysimeters

measure evaporation by monitoring the weight changes within the

enclosed soil volume while the non-weighing, 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


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.


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

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


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. Direct-weigh 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 direct-weigh 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 3-1).

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 soil-to-soil 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 soil-to-soil 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 epoxy-based 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 20-cm steel I-beams (W8x15) and extended the


entire length past the edge of the soil container (Figure 3-2). To

alleviate the need for underground access to the load cells supporting

the lysimeter, two 2-cm 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 25-cm steel I-beam

(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 I-beam. 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 10-cm

reinforced poured-in-place concrete floor was constructed for the

installation of the lysimeter at the IREP (Figure 3-2). 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


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


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


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 factory-installed 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 in-situ 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


surface). The normalized output was then recorded as the weights were


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 3-1). 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 3-3.

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-


constantan) placed throughout the soil profile. Probes were

constructed using 2-cm PVC pipe with small holes drilled radially into

the pipe at specified intervals (Table 3-2). 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

multi-pair 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


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


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

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 back-scattering 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 back-scattered 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 R-square 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


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 in-situ 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


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 3-2). 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


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 warm-up, 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


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 4-minute count was taken. In the

calibrate mode of the timer, the meter indicates the count per minute

for the 4-minute 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 one-minute 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

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 3-1).

Pwet = A B ln(CR) (3-1)


Pwet : wet bulk density of the soil [g-cm-3
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) 3-2

8 : volumetric water content [cm3cm-3]

Pwet : wet bulk density of soil [gcm-3]

Pdry : dry bulk density of soil [gcm-3]


P : density ofyater [g-cm3]
= 1.00 [g.cm ]

Combining equations (3-1) and (3-2) yields a relationship for
volumetric water content in terms of the count ratio and the dry bulk

density of the soil.

6 = A + B-ln(CR) + C-Pdry (3-3)

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 (3-1) and (3-3) 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.33-1n(CR) 0.89-Pdry (3-4)

The coefficient of variation (R2) for the regression shown in

equation (3-4) was 0.888. Measured volumetric water content when

plotted against the water content estimated by the calibration

(equation 3-4) should lie about a line with a slope of unity and an

intercept of zero (Figure 3-4). Examination of Figure (3-4) indicated

that relatively large errors could be associated with the exact
estimate of volumetric water content. Therefore, the gamma probe


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 (3-5) 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

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 3-3). 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 3-6). The error bars shown in Figure 3-6

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) PDP-11/23 mini-computer equipped with an

IEEE-488 interface board and utilizing the RT-11 operating system. The

RT-11 operating system provided several system subroutines which were

primarily for scheduling the data collection as well as other time


manipulations. The main program and associated subroutines were

written in FORTRAN-66. 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


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 PDP-11/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 CR-207 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 non-condensing 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 3-4. 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


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 ten-minute 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 day-of-year, time-of-day, the

normalized output for each load cell (mV/V), the 10-minute total net

radiometer output, the number of net radiometer readings during the

10-minute interval, the status flag, and twenty temperatures

corresponding to the depths indicated in Table 3-2.

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

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs