Illllllllllllllllilliilmlllll ,~ -- -

31293 10731 7186 f
IHESIS : Ligfififii’
g,"

This is to certify that the

thesis entitled

A TWO-DIMENSIONAL SOIL PHYSICAL MODEL
OF MOISTURE AND TEMPERATURE INTERACTIONS
IN THE SEED ZONE USING A FINITE ELEMENT FORMULATION

presented by

Don E. Holzhei

has been accepted towards fulfillment
of the requirements for

Ph. D. degree m Agricultural
Engineering

W

Major professor

 

 

 

Date August 7, 1984

 

0-7639 MSU is an Affirmative Action/Equal Opportunity Institution

 

MSU

LIBRARIES
m

‘-

 

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINE§ will
be charged if book is
returned after the date
stamped below.

 

 

mfg-.1'afigh'

 

 

 

A TWO-DIMENSIONAL SOIL PHYSICAL MODEL
OF MOISTURE AND TEMPERATURE INTERACTIONS
IN THE SEED ZONE USING A FINITE ELEMENT FORMULATION

by

Don Earl Holzhei

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements

for the degree of
DOCTOR OF PHILOSOPHY

Department of Agricultural Engineering

1984

ABSTRACT

A TWO-DIMENSIONAL SOIL PHYSICAL MODEL OF MOISTURE
AND TEMPERATURE INTERACTIONS IN THE SEED ZONE
USING A FINITE ELEMENT FORMULATION

by

Don E. Holzhei

The objective of this study was the formulation of
finite element model that could be used to analyze the mois-
ture and temperature regime around a planted seed periodi-
cally for three or more days after planting. The model
developed from this study was unique due to the combination
outlined below:

1. Studied simultaneous moisture and temperature
interactions

2. Incorporated soil air mass flow convection with
diffusion mechanisms

3. Exhibited two-dimensional and time-dependent
regimes

4. Considered heterogeneous soil grid with varying
structure and texture

5. Included soil moisture hysteresis

6. Interfaced with a diurnal meteorological model
imposed on the soil surface

7. Utilized the Galerkin approximation to the finite
element numerical solution

8. Applied to seed furrow geometry and tillage
practice studies

The grid was composed of triangular elements of varying
size, extending horizontally from the furrow centerline to a
point midway between rows, and vertically to a depth of 30
cm. The seed was considered as a passive point in the grid.
The computer program was written in FORTRAN for an IBM
360/40 computer. The model was validated using published
field data.

Simulations were run studying the following parameters:

1. Sandy loam vs. clay loam

2. Tilled vs. undisturbed soil profile

3. Bulk density comparisons in the tilled soil layer
(representing different tillage treatments)

4. Seed furrow geometry
5. Average May vs. hot June conditions for Michigan

This study concludes with suggestions regarding the
need for future work. Refinement of this model and its
expansion to include a broader scope of conditions is pre-
requisite to linking this two-dimensional model to a three-
dimensional seed-soil interface model. The latter is a
vital component of a larger seed germination model proposed
in a long-range plan.

Approved: /6/

 

 

Major Professor

 

ACKNOWLEDGMENTS

The author sincerely appreciates the guidance,
patience, and cooperation over the years of his committee
chairman, Dr. Thomas Burkhardt. His gratitude is also ex-
tended to his committee members Dr. G. Merva, Dr. L.
Segerlind, Dr. A. Smucker and Dr. A. Srivastava for their
time and counseling.

The author is indebted to Betty Nicholson and Susan
Schroeder of Computer Services at Delta College for their
tireless processing of computer runs throughout the develop-
ment. Acknowledgment is due Marnie Laurion of MSU for the
final typing of this dissertation.

Finally, the author wishes particularly to thank his
wife, Kathleen Ann, for her encouragement and support over a
period of eight years, without who's selfless devotion this
dissertation quite likely may not have been brought to
fruition. .

ii

TABLE OF

LIST OF FIGURES . . . . . . .

LIST OF TABLES . . . . . . .

LIST OF VARIABLES . . . . . .

1.

INTRODUCTION . . . . . .

1.1
1.2
1.3
1.4
1.5

Computer Modeling of

Long Range Goal . .

Objective of this Study

Soil Model . . . . .

Summary . . . . . .

CONTENTS

Page
O O O O O O O O O O O O 0 viii
O O O O O O O O O O O O 0 x
O O O O O O O O O O O O O Xii

o
o
o
o
o
o
o
o
o
o
\l

O O O O O O O O O O O O O O O 7

O O O O O O O O O O O O O O 0 10

THEORY AND LITERATURE REVIEW . . . . . . . . . . . . . 11

2.1
2.2

2.3

Moisture and Heat in

the Soil . . . . . . . . . . ll

Explanation of Philip-DeVries Equations . . . . . 12

2.2.1 Moisture Equation . . . . . . . . . . . . . 12

2.2.2 Temperature Equation . . . . . . . . . . . 13

Assumptions . . . .

O O O O O O O O O O O O O O 0 13

2.3.1 Assumptions Made by Philip and DeVries :. . 13

2.3.2 Additional Assumptions and Considerations

Simplified Equations

15

O O O O O O O O O O O O O O O 17

Theory Incorporation of Air Mass Effects . . . . . 19

2.5.1 Air Velocity Due to Pressure Gradients . . 19

2.5.2 Air Velocity Due to Density Gradients . . . 21

2.5.3 Modification of Heat Equation . . . . . . . 24

2.5.4 Modification of Moisture Equation . . . . . 27

2.5.5 Evaluation of

iii

the Velocity Gradient . . . . 30

2.6

summary O O O O O O O O O O O O O O O O O

MODEL FOR NUMERICAL SOLUTION . . . . . . . . .

3.1
3.2
303

3.4
3.5
SOIL
4.1

4.2

4.3

Background of Finite Element Modeling . .
Model Development . . . . . . . . . . . .
Integration of Terms in the Model . . . .
3.3.1 Area Integrals . . . . . . . . . .
3.3.2 Line Integrals . . . . . . . . . .
3.3.3 Time Integrals . . . . . . . . . .
Set of Algebraic Equations . . . . . . . .
Summary . . . . . . . . . . . . . . . . .
PARAMETERS . . . . . . . . . . . . .1. . .
Moisture Diffusivity De . . . . . . . . .
4.1.1 Moisture Characteristic Slope 86
4.1.2 Hydraulic Conductivity K . . . . .
4.1.3 Hysteresis EffeCts . . . . . . . .
Moisture Diffusivity D1. . . . . . . . . .
4.2.1 Factor f . . . . . . . . . . . . .
4.2.2 Moisture Contentelk . . . . . . .
4.2.3 Atmospheric Diffusion Da‘m . . . .

4.2.4 Air Pore Temperature Gradient
Ratio (VT)a/(VT) . . . . . . . . .

4.2.5 Theory Modification . . . . . . . .
Thermal Conductivity A . . . . . . . . . .
4.3.1 Soil Temperature Gradient Ratios ks

4.3.2 Temperature Gradient Ratio of the
Dispersed Fluid kav or k1 . . . .

4.3.3 Thermal Conductivity of Air-filled
Pores xa.‘, O O O O O O O O O O O O O

iv

31
32
32
34
47
47
48
49
50
56
57
57
57
58
59
60
61
61

62

62
63
64
65

66

67

4.7

Heat Capacities C and Cg .3. . . .
Vapor Content Q9 . . . . . . . .
Velocity Gradient aV/az . . . . .

4.6.1 Pressure Changes Components
3Vp/32 . o o o o o o o o

4.6.2 Density Changes Component
- BVd/Bz o o o o o o o o o

smary O O O O O O O O O O O O O

SURFACE CONDITIONS . . . . . . . . . .

5.1
5.2
5.3

Evaporation at the Surface . . . .

.Heat Flux at the Surface . . . . .

sumary O O O O O O O O O O O O O

COMPUTER IMPLEMENTATION . . . . . . . .

6.1
6.2
6.3
6.4
6.5
6.6

Finite Element Grid of Soil Furrow
Boundary and Initial Conditions .
Specific Inputs to the Model . . .
Computer Flow Diagram . . . . . .
Computer Program Listing . . . . .

Summary . . . . . . . . . . . . .

PARAMETER SENSITIVITY ANALYSIS . . . .

7.1

Profile

Soil Moisture Diffusivity De, Matric
Potential w, and Hydraulic Conductivity K

Soil Moisture Diffusivity DT . .
Soil Thermal Conductivity A . . .

Soil Heat Capacity C . . . . . .

Velocity Gradients 3Vp/az and an/az of

sail Air O O O O O O O O O O O O

Evaporation and Soil Heat Gain at the

surface O O O O O O O O O O O O O

68
69
70

70

74
75
76
76
80
82
84
84
88
89
89
90
91
92

92
97
99
102

102

103

8.

MODEL VAL I DAT I ON O O O O O O O O O O O O O O

8.1 The Bruce Study . . . . . . . . . . . .

8.2 Model Refinements . . . . . . . . . . .

8.3

8.2.1
8.2.2
8.2.3
8.2.4
8.2.5

Finite Element Formulation . . .
Soil Air Mass Flow . . . . . . .
Surface Conditions Model . . . .
Bruce Data . . . . . . . . . . .

Moisture Diffusivity Due to Vapor

Comparisons Between the Model and
PubliShEd Data O O O O O O O O O O O O O

804 summary O O O O O O O O O O O O O O O O

SIMULATIONS O O O O O O O O O O O O O O O O O

9.1 Model Parameters . . . . . . . . . . . .

9.2

9.1.1
9.1.2
9.1.3
9.1.4

Soil Profile . . . . . . . . . .
Undisturbed Soil . . . . . . . .
Tilled sail O O O O O O O O O O O

Meteorological Inputs . . . . . .

Simulation Results . . . . . . . . . . .

9.2.1
9.2.2
9.2.3
9.2.4
9.2.5

9.2.6

Temperature and Moisture in
Undisturbed Metea Sandy Loam . .

Temperature and Moisture in
Undisturbed Brookston Clay Loam .

Comparisons Between Undisturbed
and Tilled Metea Loam . . . . . .

Tilled Density and Tillage Passes
Comparisons in Metea Loam . . . .

Minimum Tillage Furrow Shape
Comparisons . . . . . . . . . . .

Varying Meteorological Conditions

9.3 Summary . . . . . . . . . . . . . . . .

vi

106
106
107
110
110
111
111

113

114
116
117
117
117
117
122
124
127

127

130

133

136

136
142
142

10. CONCLUSIONS

10.1 Summary

AND FUTURE WORK . . . . . . . .

10.2 Conclusions . . . . . . . . . . . . . .

10 O 3 Future work O O O O O O O O O O O O O O

10.3.1
10.3.2
10.3.3
10.3.4
10.3.5
10.3.6

Second Validation . . . . . . .
Refinement of Parameters . . . .
Model Expansion . . . . . . . .
The Three-Dimensional Seed Model
The Three-Dimensional Soil Model

The Seed-Soil System Model . . .

LIST or REFERENCES . . . . . . . . . . . . . . .

APPENDIX A-COMPUTER FLOW DIAGRAM

APPENDIX B-COMPUTER PROGRAM LISTING

APPENDIX C-CONVERSION TABLE

APPENDIX D-METEOROLOGICAL INPUT CALCULATIONS

vii

145
145
146
147
147
148
149
150
150
151

152

Figure
101

1.2
2.1
2.2
3.1
6.1
6.2
6.3
7.1

8.1
8.2
8.3
9.1
9.2
9.3

9.4

LIST OF FIGURES

Seed-Soil System Model for Germination . .
Soil Model Decision Tree Diagram . . . . .
Infinitesimal Element of Soil for Heat Flow
Infinitesimal Element of Soil for Mass Flow
Gray and Pinder Model . . . . . . . . . .
Finite Element Domain of Soil Profile . .
Finite Element Grid for Validation . . . .
Finite Element Grid for Simulation . . . .

Moisture Characteristic Curves for
Adelanto Loam . . . . . . . . . . . . . .

Hydraulic Conductivity for Adelanto Loam .
Moisture Diffusivity Db for Adelanto Loam
Thermal Conductivity of Quartz Sand at 20°C
Thermal Conductivity of Quartz Sand at 40°C

Daily Evaporation for a Loam in Binghamton,
NY, in June . . . . . . . . . . . . . . .

June 18 Temperatures comparing Model with
Bruce Data O O O O O O O O O O O O O O O O

June 18 Moistures comparing Model with
Bruce Data O O O O O O O O O O O O O O O O

June 18 Noon Moistures comparing Liquid and
Vapor Diffusivities . . . . . . . . . . .

Moisture Characteristic Curve for Metea
sandy Loam O O O O O O O O O O O O O O O O

Moisture Characteristic Curve for Brookston
Clay Loam O O O O O O O O O O O O O O O O

Effect of Soil Structure on Moisture
Characteristics . . . . . . . . . . . . .

Ghosh Argument for One ¢b . . . . . . . .

viii

24
28
35
85
86
87

93
94
95
100
101

104

108

109

115

122

123

9.5

9.6

9.7

908

9.9

9.10
9011

9.12

9.13
9.14

SO15

9.16

A.l

D.1

Temperatures in Undisturbed Metea Sandy
Loam for May . . . . . . . . . . . . . . . .

Noon Moistures in Undisturbed Metea Sandy
Loam for May O O O O O O O O O O O O O O O O

Temperatures in Undisturbed Brookston Clay
Loam for May O O O O O O O O O O O O O O O O

Noon Moistures in Undisturbed Brookston Clay
Loam for May O O O O O O O O O O O O O O O O

Moistures at Selected Depths in Metea Sandy
Loam comparing Undisturbed and Tilled . . .

Noon Moisture Profiles in Metea Sandy Loam .

Third Day Temperatures in Metea Sandy Loam
comparing Undisturbed and Tilled . . . . . .

Third Day Moistures in Metea Sandy Loam
comparing Tilled Densities . . . . . . . . .

Minimum Tillage Furrows Modeled . . . . . .

Moistures at Seed Point in Metea Sandy Loam
comparing Furrow Shapes . . . . . . . . . .

Third Day Moistures in Metea Sandy Loam
comparing Meteorological Conditions . . . .

Third Day Moistures in Metea Sandy Loam
comparing Meteorological Conditions . . . .

Computer Flow Diagram

Global Irradiance Cycle Computation

ix

128

129

131

132

134

135

137

138
139

140

143

144

Table
4.1

4.2

9.2

9.3

9.4

905

LIST OF TABLES

" .‘o. .

Moisture contents at K a Ks/lOOO for several

soils (Mualem (1978)) . . . . . . . . . . . .

Thermal conductivities of several elements,
J/m .Sec . 0c O O O O O O O O O O O O O O O O 0

Soil temperature gradient ratios . . . . . . .
Heat capacities, J/m3°°C . . . . . . . . .

Air permeabilities, K*, m/sec . . . . . . . .
Volume surface mean diameters, mm . . . . . .

Air pressure wave measurements (Farrell, et
alO (1966)) O O O O O O O O O O O O O O‘O O O

Surface roughness, 20, cm . . . . . . . . . .
Soil surface albedo . . . . . . . . . . . . .
Depth of influence . . . . . . . . . . . . . .
Boundary and initial conditons for the model .
Soil moisture properties . . . . . . . . . . .
Moisture diffusivities,inZ-secjl/1°C—or-6). .

Velocity gradients, m/sec-m . . . . . . . . .
Aerodynamic resistance for AT 8 4°C . . . . .

Soil conditions in Northern locations
pertinent to model . . . . . . . . . . . . . .

Meteorological statistics from Bishop
Airport (1941-1970) 0 o o o o o o o o o o o o

Diurnal meteorological statistics for
1981 from Bishop Airport . . . . . . . . . . .

Mean weekly solar radiation over East Lansing,
langley/day . . . . . . . . . . . . . . . . .

Meteorological inputs used in the simulations

Page

62

66
66
69
71
72

74
79
81
84
88
97
98
103

105

121

125

125

126

127

B.1

C.1

Computer program

Conversion table

xi

al

3:119

atm

aw/ae

LIST OF VARIABLES
volumetric air content, In3 of air/m3 total volume
albedo, reflectivity coefficient of soil surface
amplitude of pressure waves of air moving

horizontally

specific heat of liquid water, J kg.1 C-l

specific heat of water vapor at constant pressure,

J kg.l C'-1

volumetric heat capacity of moist porous medium,

J m—3 C"1

volumetric heat capacity of soil, liquid, and gas
constituents of a soil medium

molecular diffusion coefficient of water vapor in

air, 4.42x10-4T2°3/P, cm2 sec-1, or 5.801x10-11T2°3,

m2 sec-1

thermal moisture diffusivity, DTl+DTv' 1112 sec".1 C-l

liquid diffusivity due to temperature gradients,

2 -1

kyu, m sec-l C

vapor diffusivity due to temperature gradients,

fDatmv8h(VT)a/p1VT, m2 sec.1 C-1

2 -1

isothermal moisture diffusivity, D61 + D , m sec

6V

liquid diffusivity due to moisture gradients, Raw/39

m2 sec-1

1'

vapor diffusivity due to moisture gradients,

2 -l
aaDatmvgpv(aw/ael)/leT, m sec

lepe of soil moisture characteristic curve, m H20

evaporation rate, msec 1

xii

mass evaporation rate, kg m-Zsec-l

factor accounting for liquid island effect,
describing effective space available for flow,
exclusive of tortuosity,

S for el<e1k

a+ael/(S-elk) for 6 >6

1 1k

acceleration due to gravity, 9.81, msec"2

relative humidity

absolute humidity of air, kg m-3

saturation absolute humidity at the surface

temperature, kg m—3

absolute humidity of air at the soil surface,

kgm-3

mechanical equivalent of heat, 4.18x107, erg cal-

unit vector in vertical direction

unsaturated hydraulic conductivity, m sec-1

hydraulic conductivity at saturation, m sec"-1

0 cm-1

6

air permeability, cm sec-l/cm H2

heat of vaporization, 2.45115x10 l

.Jkg-
partial pressure of water vapor, mm Hg

total gas pressure, mm Hg

sensible heat flux to the air, W In-2

total global short wave irradiance from sun and

atmosphere, W m-2

incoming long wave radiation from sky, W m-2

net irradiance, W m 2

sensible heat flux in soil, W m-2

xiii

l

6 1 1

gas constant of water vapor, 4.615x10 erg g- C-

or 0.11 cal gfllc-l

adiabatic or neutral value of Rc’ sec m-1

aerodynamic resistance, sec m-l
Richardson Number

porosity, m3/m3
wind speed, measured at 2m elevation, m sec-l
stability correction factor

time, sec

temperature, °K or °C

air temperature, °C

surface temperature, °C

vertical air mass velocity due to density gradients
vertical air mass velocity due to pressure gradients
horizontal coordinate, m

vertical coordinate, m

surface roughness, m

tortuosity factor for diffusion of gases in soils,

0.67

3 3 -1

dpo/dT, 1.05x10’ @ 20c, kg m' c
or = dimensionless empirical constant in soil

w-e relationship

temperature coefficient of surface tension of water,
-2.09x10'3 @ 20c, c"1
emissivity of soil surface (dimensionless)
total volumetric content of moisture, 61 + 6v,
m3/m3 soil volume

xiv

1,2,3

volumetric moisture content at air entry point (at
"knee" of w-e curve)
volumetric liquid content, m3/m3 soil volume

value of e at which "liquid continuity" fails

1
moisture content at a point between wilting and

hygrosc0pic levels where soil is considered dried
volumetric vapor content, m3 of precipitable water/m3
soil volume

thermal conductivity of soil, weighted average over
its constituents, inkiAi/ZXiki, (where xi = volume

fraction, k1 = ratio of temperature gradient across

ith constituent to overall gradient), J m”1 sec-1

C-l

. . . -1 -1 -1
thermal conduct1v1ty of dry air, J m sec C
thermal conductivity of air-vapor mixture,

-1 -l -l
Aa+lv, J m sec C

thermal conductivity of liquid, J mnlsec-J'C-l

1 1 -1

thermal conductivity of soil, J m- sec- C

thermal conductivity of vapor, apparent in an air
filled pore due to vapor diffusion, hvBLD

J m-l'sec"]'c-l

atm'

% sand, silt, and clay, respectively, in soil
volume

P/(P-p), mass flow factor due to diffusion within
an air pore‘

soil bulk density, kg m-3

density of liquid water, kg mm3

XV

3

p0 = density of saturated water vapor, kg m-

pv = density of water vapor, poexp(ug/RT), poh, kg m-3

o = Stefan-Boltzmann constant, 5.67x10-8, W m-zK-4

w = soil moisture potential, matric, m H20 bars or
atmospheres

we = air entry value of matric potential, obtained from
w-e curve, evaluated at 6e

wsv = matric potential at the soil surface, m H20

w = frequency of air pressure waves moving horizontally,
sec-l

(VT)a = average temperature gradient in air-filled pores,
C m-1

xVi

1. INTRODUCTION

1.1 Computer Modeling gf Soils

The computer has become an extremely powerful tool in
design and research. It allows the engineer to analyze a
part or product providing information which may not have
been obtainable before, or only obtainable with a costly
prototype exposed to extensive testing. Software programs
using numerical methods such as finite elements are used to
perform stress analysis, vibration analysis, or kinematic
analysis. Manufactured products and natural processes alike
have been modeled to predict their response to an external
stimulus through the method of simulation. This adds to the
understanding of complex interrelationships and provides a
better basis for decision-making for equipment or process
managers.

One area which can benefit substantially from computer
modeling is that of soils research. Specifically, a better
understanding of the immediate environment around a seed
after it has been planted and during the germination and
emergence periods could:

I. serve as an educational tool for study by students
from many disciplines,

2. aid the researcher in the study of different tillage
methods such as minimum-till, ridge tilling or slot mulch
tilling,

3. aid the farmer in crop residue management,

2

4. aid the machinery manufacturer in tillage equipment
design.

Agricultural and engineering students can all benefit
from such a soil model by observing moisture and temperature
changes in the soil profile which could not be fully appre-
ciated in a lecture on the subject nor by extensive soil
testing in the lab.

Researchers should be able to determine, for a given
climate, seed and soil combination, whether the seed dis-
plays more or less viability during the germination process
for one tillage treatment over another. Studies could be
made on the immediate seed environment, how it is affected
by the tillage treatment and climate over a period of days,
and how it in turn affects the seed and its chances for
survival. Complex soil moisture-temperature interactions
can be studied, such as the effect of a tilled soil layer on
the impedance to moisture flow as the layer heats vup and
dries, the effects of bulk density on heat and moisture, and
hysteretic effects due to diurnal variations.

The farmer could use the model on a home computer for
managerial decisions on his operation. Example could be "If
I plant corn shallow with conventional tillage, expecting a
rain, but instead it stays dry for six days, what happens to
the moisture at the seed level, and will I still get an
acceptable germination?" ”If I plant sugar beets now, and
it freezes two nights in a row, what can I expect for a

germination response?" "With the present soil

3
condition and weather, how long can I delay planting and
still be assured of an acceptable stand?"

The planter/tillage equipment manufacturer could use
the model to study furrow depths and shapes, soil aggrega-
tion, and seed placement relative to the furrow to gain
insight as to which tillage implement combinations are best
suited for specific regions, soils, and crop cultures. With
the model the manufacturer could determine whether furrow
shape or width is the prime factor in moisture retention at
the seed level and what bulk density above and below the
seed is optimum for seed germination.

Much has been published on soils, and on crop response
to soil, meteorological and tillage conditions. Through the
decades the body of knowledge of soils and crops has grown
considerably, and recently some models of systems and sub-
systems have been proposed. Some of these models have been
parochial in nature, concentrating on a subsystem artifi-
cially isolated from other pertinent factors in the total
system. Hillel (1977) made detailed computer simulations of
heat and moisture interactions in the soil, but confined
his study to one-dimensional flow using relatively large
soil increments in the model. More recently Cruse, et a1.
(1980) develOped a model to predict temperatures at the 5-
cm. depth as affected by tillage, but that model did not

incorporate direct interactions with moisture. Gray and

Pinder (1974) applied the finite element analysis to tran-
sient groundwater flow, but it was confined to isothermal

conditions.

Other models have been broadly-based crop growth models
to predict growth response to macro-meteorological and soil
conditions, but do not incorporate the intricate heat-mois-
ture-oxygen-crop interactions in the soil. Examples of
these are the CORNGRO (Tscheschke and Gilley (1979)) and

SORGF (Arkin, et a1. (1980)) models.

L2 Lora Ram: 9211.

The long range goal surrounding this study was to
develop a model of the soil environment around a planted
seed with sufficient detail to understand all pertinent
climate-soil-seed interactions in the system. The intended
use of the model would be to study the seed-soil interac-
tions for the duration of the germination period to predict
seed viability in response to factors such as:

1. soil texture

2. meteorological conditions

3. tillage methods and soil structure

4. cropping management

Thus this model could be used to determine whether a
given tillage practice, under specified weather conditions,
would provide a suitable environment for seed germination.

The emergence of the seedling is also crucial to

developing a healthy crop stand. This emergence period

would require a seedling model different from the seed model
and would necessitate incorporation of soil chemistry for
surface crusting conditions.‘ Since germination and emer—
gence mechanisms are significantly different and do not
occur at the same time, the author decided a completely
separate emergence system model could sometime in the future
be interfaced with the germination system model of Figure
1.1.

Figure 1.1 represents the long range goal. A two-
dimensional soil physical model views the soil profile in a
vertical plane sliced across the furrow with the seed con-
sidered as a passive point in the profile. Meteorological
diurnal data from the meteorological model are superimposed
onto the soil surface through a soil surface interface
model. At each time step, physical variables such as mois-
ture and temperature are computed throughout the profile.
The computed variables near the seed are then transferred
into the grid of a 3-D soil model which forms a shell around
an active spherical seed model. The 3-D soil model inter-
acts with the seed model within each time step. It extends
only far enough to incorporate the seed's sphere of influ-
ence, and includes a seed-soil interface boundary across
which moisture, heat and oxygen exchange. The 3-D soil
model thus models only the immediate environment around the
seed, whereas the 2-D soil model models the entire soil

profile.

cowuocwfiuoo new Hone: Eoum>m aflomlpoom .H.H ousmflh

 

 

 

 

00mm — 0 0“
Home Hwom cw xomn 9
a umsnouiom _ uuo>coo

 

 

    
 
    

 

   

Hobo: mEOADHUcoo
aim aflomipoom mumpcson
aw anewuoououcH can

  
    

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ouzmeoo noun Hufiuwcw 3oz
med» Hum
mooa use
a - - noun
Mono: ouonmm o» comm Home mafia co>wm now ~o Hobo:
anon cum 2 uuw>coo .9 .o ousaeoo Haoamssa Hwom aim
Honor , A A x
wmmnw.m.._
aim boom 1 x 9

 

deco: ,
HmoflmoHouoouoz

 

7

After the 3-D soil model has equilibriated for the
given time step, the new physical variables are re-intro-
duced back into the 2-D soil model as initial and boundary

conditions for the next time step.

1.3 Objective 2£ This Study

The system model in Figure 1.1 represents many man-
years of development. To fulfill the requirements of a
Ph.D. dissertation within a reasonable time span, it was
necessary to concentrate on developing a portion of the
system model. This portion would require the following:

1. It must be self-sustaining, that is, embody a
sufficiently comprehensive system within itself to be veri-
fiable with actual field data.

2. It must be easily interfaced with the rest of the
total system as other models are developed.

The objective of this study was to develop, verify and
experiment with the 2-D soil model coupled with the
meteorological and surface interface models, considering
the seed as a passive point in the profile. The scope of

this model development is clarified in the next section.

112 Soil Model

In order to narrow the scope of development of the 2-D
soil model to a manageable level, certain assumptions and
simplifications had to be made. Figure 1.2 delineates these
decisions in a decision tree diagram form. Remarks con-

cerning the decisions are numbered along the side of Figure

 

SOIL MODEL

 

 

 

 

SURFACE PLANTING DEPTH Ap HORIZON 30-40 cm
(5 cm) (18-20 cm) DEPTH
IIJP
an s-o

 

BIOLOGY PHYSICS CHEMISTRY

MOISTURE

 

     

 

TEMPERATURE

[OXYGENI ISOIL IMPEDANCE]

 

 

VAPOR LIQUID lCONDUCTIONJ [CONVECTION]
FLOW FLOW

 

 

 

 

IIMOISTURE-TEMPERATURE INTERACTIONS I

 

 

 

 

 

 

 

 

J7 l
IRREVERSIBLE MECHANISTIC-
THERMODYNAMICS DARCY/CONTINUITY/
(Taylor & Cary) CONSERV. OF ENERGY

EQUATIONS
(Philip & DeVries)
v

 

ISOTROPIC ANISOTROPIC
HOMOGENEOUS HETEROGENEOUS

STEADY STATE

 

 

 
   

HYSTERESIS NON-
HYSTERETIC

Figure 1.2. Soil Model Decision Tree Diagram

 

N
o

1.2 and listed below.

1. The entire soil profile must be considered to a
depth well below the seed where temperature and moisture
gradients are no longer significant. The seed environment

will be influenced not only by soil conditions above it, but

also by those below such as plow pan or B1 horizon texture
significantly different from the AP horizon.

2. Soil chemistry and microbiology definitely play a
role in seed-soil interactions during germination. But
processes in the 2-D soil model alone are more apt to be
dominated by physical factors, thus the model was confined
to a physical system.

3. Moisture and temperature are Chosen as the two
physical factors to be studied together. There has been
considerable research done on these factors over the years,
and they often dominate other factors in the sOil. Oxygen
can be added when considering seed-soil interactions in the
3-D model, and soil impedance would come into the separate
emergence system model.

4. The mechanistic phenomena the Philip/DeVries equa-
tions are based on are more easily understood and have con-
siderably more research to back them up than do the
Taylor/Cary equations.

5. Isotropic conditions can be assumed throughout the
soil profile except at the interface between tilled and
undisturbed soil. Differences in soil texture and structure

throughout the profile makes it heterogeneous.

10

1.5 Summary
This chapter outlined a long range goal for seed ger-

mination modeling, and the short range objective of this
- dissertation to develop the 2-D soil model. Chapter 2 lays

the theoretical foundation for this soil model. Chapter 3

develops the partial differential equations into a matrix
of algebraic equations suitable for modeling on a computer.
The soil model calls for soil parameters. Expressions for
these are developed in Chapter 4 to the point of requiring
model inputs which are easily measured. The meteorological
model and soil surface interface model with surface
evaporation and heat flux are developed in Chapter 5. Grids.
for the finite element solution and initial and boundary
conditions are outlined in Chapter 6 for computer implemen-
tation. Chapter 7 discusses the sensitivity of soil para-
meters to computer inputs and surface conditions. Chapter 8
compares the model results with field data for verification
of the model. Chapter 9 displays a variety of simulation
results demonstrating the versatility of the soil model.
Chapter 10 presents conclusions and proposes future work
beyond the scope of this dissertation which will meet the

long range goal outlined in Figure 1.1.

2. THEORY AND LITERATURE REVIEW

The theory of moisture and heat interactions in the
soil was first advanced in the 1950's by Philip and DeVries
(1957) and DeVries (1958, 1963). They made use of the con-
tinuity equation, Darcy's law and the diffusion equation
applied to separate liquid and vapor mechanisms to obtain
two diffusion-type equations as follows:

Moisture:

(1 + Dev/(avDa ) - pV/ol)861/3t + ((S-el)hB/ol)3T/3t = [1]

tm

V(Deve1) + V(DTVT) + aK/az

Temperature 3

.-1
(c + L(S-61)h8)8T/3t + (LolDev/(avDatm) - LoV + 013 g

(u-T3m/3T)) sol/at = V(AVT) + LolV(DeVV61) + plclupelvel

+ DTIVT + K§)vm) + plcp((DeVV61 + DTVVT)VT) [2]
where D6" = vapor diffusivity due to moisture gradients
a = tortuosity factor for diffusion of gases in soils
v = mass flow factor due to diffusion within an air
pore
Datm = molecular diffusion coefficient of water vapor
in air
pv = density of water vapor
of = density of liquid water
61 = volumetric liquid content
h = relative humidity
S = porosity
D = isothermal moisture diffusivity

11

12

DT = thermal moisture diffusivity

K = unsaturated hydraulic conductivity

2 a vertical coordinate in soil profile

C = volumetric heat capacity of moist porous medium

L = heat of vaporization

B = change in density of saturated water vapor with
temperature

w = soil moisture matric potential

c1 = specific heat of liquid water

c = specific heat of water vapor at constant

p pressure

1 = soil thermal conductivity

2.2 Explanation g; Philip-DeVries Equations

2.2.1 Moisture Equation.

86
L
The terms (Dev/ownatm - ov/ozygg-
and ((s-e£)hB/o£) 3%

of equation [1] are products of differentiation to satisfy
the continuity principle and represent the time rate of
change of the volumetric vapor content 6v. This is appli-
cable where Changes in liquid and vapor content are of the
same order (dry soils, high temperatures). The remaining
term on the left,ae£/at, represents the time rate of change
of volumetric liquid content of the differential element,
also satisfying the law of conservation of mass.

The first term on the right, V(De V9£)r18 an expression
for the total moisture flux due to liquid moisture gradients
alone. The second term, V(DTVT), expresses the total mois-
ture flux due to thermal gradients alone. Both fluxes are

assumed additive, and are based on diffusivities which are

13

in turn based on the assumption that the liquid and vapor
mechanisms within each gradient influence are additive. The
last term on the right, aK/az, evolves from a consideration
. of the gravity effect on moisture movement in the expression

for total moisture potential, ¢ - w matric - 2.

2.2.2 Temperature Equation. This equation is based on the
conservation of heat principle and recognizes the total heat
content of a volume of soil is made up of the heat capacity
of the soil particles plus the latent heat of the water
vapor plus the sensible heat of the liquid and vapor separ-
ately minus the differential heat of wetting.

The terms 148-e£)h8%%
ae
-Lp, __£_
tm v at
of equation [2] represent the heat transfer due to moisture

and (LpiDev/avna

Changes between the liquid and vapor phases. The term
(pzj'lg(w-Taw/3T));%£: represents the differential heat of
wetting. The term V(AVT) represents heat transfer due to
pure conduction. The‘ term Lo£V(DevV6£) represents the
transfer of latent heat by vapor movement. The terms

°z°£(‘°ezV°z + DTL VT + KK)VT)andp£Cp((DevV6£ + DTVVT)VT)

represent the transfer of sensible heat in the liquid and

vapor form, respectively.

2.3 Assumptions
2.3.1 Assumptions Made by Philip and DeVries.

a. Moisture vapor flows through dry soil by a series-

parallel path around and through liquid islands via

14

an evaporation-diffusion-condensation-evaporation
mechanism, accounted for by the factor f in the

expression for the vapor diffusivity

DTV = fDatm vBh(VT)a/O£VT.

Liquid and vapor diffusivities due to a potential

gradient are additive.

Moisture flux due to a temperature gradient is
additive with that due to a moisture gradient.
Temperature gradient across an individual air pore
may be appreciably higher than the average across
the soil matrix at that point. This is
incorporated as a pore-soil temperature gradient
ratio (VT)a/VT into the expression for the thermal‘
vapor diffusivity DTV.

Parameters can be evaluated and are assumed reason-
ably constant in the 10°C to 30°C range.

In the thermal vapor diffusivity term, as e in-
creases above the critical level 61k where liquid
continuity results, the effective cross-section for
combined liquid-vapor transfer will decrease
steadily.

All processes of heat transfer and sources and
sinks are evenly distributed throughout the soil
volume.

Heat transfer by convection and radiation is negli-
gible.

Heat capacity of the soil is the weighted average

15

of the capacities of the solid, liquid and gas
constituents.

Thermal conductivity of the soil is the weighted
average, accounting for particle shape, of the
several constituents. Further, the thermal conduc-
tivity of the air constituent equals that of dry
air plus that of the vapor. The vapor term takes
into account the vapor distillation due to tempera-
ture gradients in the air-filled pores. This leads
to an apparent increase in the thermal conductivity
of the air-filled pores.

Generation of heat due to viscosity of the moving
liquid is negligible.

A unique relation exists between w and e, (no

hysteresis).

2.3.2 Additional Assumptions and Considerations

There are additional simplifying assumptions recognized

by DeVries (1958) which hold true for most soil conditions

at virtually no loss to generality.

a.

Sensible heat transfer by vapor diffusion is
usually negligible. This is the last term in equa—
tion [2].

Changes of moisture contained in the vapor phase
over time (aev/at) are often assumed negligible.
This is because vapor content is usually small

compared with liquid content, and the relative

l6

humidity is essentially constant at 100% throughout
the moisture range until the wilting point is
reached. This is then true in moist soils under
moderate or low temperatures. This affects many of

the terms on the left hand side of both equations.

With the environmental and time conditions outlined in

Chapter One, further assumptions were made by the author:

Heat of wetting, significant when irrigating hot
and dry soils, is negligible under these
conditions.

Soil is sufficiently moist throughout the period so
that vapor content is very small compared with
liquid content. Thus assumption (b) above is
further justified.

Sensible heat transfer by liquid movement is negli—
gible. DeVries already alluded to this to be the
case for most conditions. It will be further jus-
tified if the moisture content is always considered
below field capacity. In this region the heat flux
due to pure conduction would be the dominant term.
Gravity effects on liquid flow are negligible. At
moisture contents below the field capacity, as the
soil dries the matric potential becomes large com-
pared to gravity potential. Even in a fairly wet
soil, consideration of only a few centimeters of
soil depth in the seed zone would produce negli-

gible hydraulic head.

17

e. Vapor diffusivity due to moisture gradients (Dev)
is negligible. DeVries (1958) showed this to be
significant only at the extreme low end of the
total moisture range, but virtually zero compared

to De throughout the useful part of the range.

1

2.4 Simplified Equations
Utilizing all of the above assumptions, the heat-moisture

equations reduce to:

Moisture: ae/at V-(DTVT) + v-(bevel) [3]

Temperature: (:3T/3t = v-(IVT) [4]

A number of investigators have scrutinized the thermal
conductivity expression proposed by Philip and DeVries
(Hadas, 1969; Kimball, et al., 1976; Hadas, 1977a, b;
Sepaskhah and Boersma, 1979). Others have tested the theory
under varying conditions and proposed further modifications
(WOodside and Kuzmak, 1958; Hadas, 1968; Laroussi, et al.,
1975; Malik, et al., 1979; Jury and Letey, 1979). In most
cases the investigations have supported the theory and have
clarified the limits or proposed modifications to diffusivi-
ties or the conductivity term.

However, Hadas (1968, 1969, 1977a) has questioned spe-
cifically the assumptions (h) and (1) listed in Section
2.3.1. The two diffusion-type equations do not predict
adequately the moisture-heat interactions under diurnal

temperature conditions in natural soils.

18

The equations are based on negligible heat convection.
However, near the surface, heat and vapor flow due to air
mass movement through the soil may be appreciable. In fact,
Ojeniyi and Dexter (1979) state that much of the drying of
tilled soil in the field occurs as a result of convective
transport of air through large pores produced by the till-
age. This air mass movement would be induced by air pres-
sure gradients arising from wind gustiness at the surface
and by air density gradients arising from reversed tempera-
ture gradients in the top few centimeters of the soil in the
evening and early morning hours.

Hadas (1969, 1977a) and others have proposed a mass
enhancement factor: fitted into the vapor thermal conductiv-
ity term to account for this: lav = Aa + 5 Av

More recently Hadas (1977a) recognized the probable
inadequacy of this method and suggested in its place a third
equation be added to the theory to account for these diurnal
effects near the surface. Similarly, Farrell et al. (1966)
had proposed a mass velocity squared term be added to the
vapor diffusivity to account for increased vapor movement
near the surface.

Both of the above methods involve correction terms to
bring theory in line with experimental data. An alternative
method is to return to the basic study of heat and mass flow
through a differential volume of soil. Below is presented
the development of this alternate method to incorporate air

mass movement effects in the surface soil layers.

19

2.5 Theory Incorporation gf Air Mass Effects

Luikov (1966) developed a set of three equations using
the thermodynamic approach to explain heat and mass transfer
through capillary-porous bodies under the influence of tem-

perature, mass and gas pressure gradients:

aT/Bt = 1(1va

_ 2 2 2 [5]
ae/at — sz e + K3V T + K4V p
Bp/Bt = K5V2p

Where the Hi terms are phenomenological transfer coeffi-
cients. These equations still do not include density gradi-
ent effects. To incorporate density effects, a fourth equa-

tion would have to be added of the type: apg/at = K V20

6
where 09,- density of the gas. 9

Since this would produce a number of transfer coeffi-
cients for which there is little if any experimental data, a
different direction is proposed.

Both pressure and density gradients give rise to verti-
cal momentum of the air-vapor (gas) mixture in the soil,
especially in tilled, porous, drier top layers of soil. If
the net velocity of the gas due to these gradients can be
calculated at any point in the soil at any time, then this
velocity can be incorporated into an expanded expression of

the original two heat-mass equations, as suggested by

several investigators.

2.5.1 Air Velocity Due To Pressure Gradients. Several

investigators have studied this phenomenon. Fukuda (1955)

20

developed a mathematical equation, based on the Darcy law
and supported by data, for calculating the transmission of
air pressure and air mass movement through the soil due to
surface pressure disturbances. His governing equation was
of the form:
splat = (k/a)a(pap/az)/az

where p = air pressure, k a air permeability, and a = volu-
metric air content, or air porosity. His direct solution
was for a vertical periodic pressure disturbance at the
surface.

Hanks and Woodruff (1958) discussed, but did not prove,
that vapor transfer would be also enhanced by a ”mixing
effect" caused by pressure fluctuations in addition to air
mass actually leaving the soil as in Fukuda's model.

Farrell, et al. (1966) expanded on Fukuda's model,
incorporating a second spatial dimension, and allowing for
variable soil depth to an impermeable layer and variable
wave length of the surface disturbance. Horizontal effects
were found to be negligible. Their macroscopic vertical air
velocity expression VEJwas:

** * 'k *
Vp = -(K A ao/e)exp(-l z)[cos(wt-y z-2ny/L )
[6]

* 4
+ r sin(mt-Y z-2wy/L )]

* * *
2 v 2 + (Zr/L )2

* 'k
where l and y are derived from: A

e
A

*
*
s/2K poy
and:

e = soil porosity

21

p - mean atmospheric pressure, cm.Hé0

L - wavelength of horizontal pressure waves, cm

y - horizontal coordinate, cm

2 - depth, cm

r - Y*/X*

K a air permeability, cm sec-llcm ago cm.-

a - amplitude of pressure waves of air moving horizontally, cm

m . frequency of air pressure waves moving horizontally, sec

t - time, sec

2.5.2 Air Velocity Due To Density Gradients. Luikov (1966)

and Hollman (1976) gave a brief discussion of free convec-

tion in enclosed spaces. Their treatment was of two hori-

zontal plates Closed around the edges with the bottom plate

heated. Earlier investigators had shown fluid flow in this!
free convection mode to be a pattern of mushroom Clouds, or

'Benard cells“, positioned perpendicular to the plates.

Fluid moves up from the hot plate inside each cell and as it

cools at the cold plate it falls back down around the out-

side of each cell. Although this continuous Circulation of

fluid "pumps“ heat across the plates by free convection,

Hollman (1976) outlined that it is often treated as conduc-‘
tion. An apparent thermal conductivity ke representing

convection was identified in the Nusselt Number Nu - ke/k,

where Nu is a function of both the Grashof and Prandtl

Numbers. Knowing geometrical and fluid parameters, this

apparent conductivity can be calculated.

Hadas (1977a) proposed that free convection in air-

22

filled. pores in soil may be similar to this process. A
Closer analysis would show that, although this may be the
dominant mechanism in wet or deep layered soils where air
. pores are isolated from each other, it is probably not the
dominant mechanism in drier, porous, well-tilled soil layers
near the surface. Here the pores would provide a more cont-
inuous medium for air mass flow. This is evidenced by field
measurements of appreciable moisture movement back towards
the surface during temperature reversals in the drier sur-
face layers.

Assuming this latter condition, density gradients would
then give rise to buoyant forces pushing the air mass verti-
cally through the soil. The following analysis is based on
this assumption.

According to Gebhart (1961), the difference between the
weight of a small volume of fluid and the buoyant force on

it due to a density gradient is the net force

F a 98*ov(T-Tw)

where: F . force/unit volume, N m“3
g - gravitational constant, m sec"2
*
B a coefficient of thermal expansion, °H 1

p - fluid density, kg m"3

T - temperature of the fluid volume under
consideration, °R

T - temperature of the total fluid medium, °H

23

The actual path of this small volume of gas follows a
tortuous route on its way to the surface. If we assume the
macroscopic kinetic and potential energy of the volume is
conserved, then a relation between the net force and
velocity can be obtained.

The change in potential energy of the gas volume across
a vertical increment of soil of depth Az can be expressed
as:

A P.E. = 98*vaTAZ
The change in kinetic energy of this same volume across

this same soil increment can be expressed as:

__ 1 2 _ 2

change in temperature across the soil increment

where: AT

Vdu a macroscopic velocity at upper layer of soil

increment
le . macroscopic velocity at lower layer of soil
increment

Equation the two energy terms and rearranging:

 

7

vdu = 1/2gB*ATAz + Val2 [71
Equation [7] requires the knowledge of the macroscopic
velocity of air at the lower boundary of the soil increment
and the temperature gradient across the increment. Thus,
the velocity of air due to density gradients will depend on
vwhere the element in question is positioned vertically with-

in the soil profile.

24

2.5.3 Mbdification of Heat Equation. The total air
velocity is, from equations (6) and (7),

_ V = Vp + ivd [8]

1 if temperature gradients are reversed (hot sub-

where i = layer)

0 if temperature gradients are normal (hot surface)

 

This is now incorporated into the heat equation by returning
to a basic study of an infinitesimal element of soil com-
posed of solid, liquid and gas (air-vapor) constituents, as

shown in Figure 2.1.

oz +an I azdz
qy
_ t /’
| 5
dx

+aqylay~dy

 

dz

q ->1 ——>qx+aqx/ax-dx

 

 

 

 

1

C12

Figure 2.1. Infinitesimal Element of Soil for Heat Flow

25

Referring to Figure 2.1., the total heat flux per unit area
of this element in direction n is:
qn a qn conduction + qn free convection + qn forced
convection
s qn conduction + qn convection
and generally, qn conduction = -AnaT/an
If the total convection term (qn convection) is abbreviated
gen, then the difference in rate of heat flow per unit area
along the x axis is:
- 3% (-Ax %% + qcx) dx

The net rate of gain in this element of volume dxdydz due to

this flux is:

3 8T
3;; ”x I; + qcx) dx ‘dydzl
3 3T 3
or 3x (Ax 8x) + 3; (qcx))dxdydz

Similar expressions can be obtained for heat flux along the

y and z axes. A general expression would be of the form
(V-(AVT) + ch)dxdydz

If only the x and 2 directions are considered, and convec-

tive currents are confined to the z direction, then the net

rate of gain becomes:

8T

3T 3 aqcz
ax) + (AZ 3;) + -—§;) dxdydz

8
(a; ” 32

X
The time rate of Change of stored energy in the element is:
m-CvDT/Dt
where m - mass of element
c - mass heat capacity of element

DT /Dt - particle derivative of temperature of element

= aT/at + uaT/ax + vaT/ay + waT/Bz

26

u,v,w, - velocity components of the element in the x,y,z

directions, respectively
The element's solid and liquid constituents are at zero
- velocity, but the gas is not. Thus, the stored energy term

can be expanded to reflect both temperature and momentum

changes. DT DT DT
S S Dt 1 1 Dt g g Dt

where s,l,g - solid, liquid and gas subscripts.
Letting p - m/dxdydz and assuming 09- constant at
these low velocities, the time rate of change of stored

energy becomes: Tl

D
olcl BE. + 0999 Dt ) dxdydz

L3

(0 c +

__E
s s Dt
and since the only velocity consideration is that of the gas.
in the vertical (z) direction, and assuming all constituents
are in thermal equilibrium at any one point, the expression

reduces to: 3T

(0 c ‘—— +

8T
8 s at Ole + ogcg (3?

é‘filf‘a’

3T
1 + w 35))dxdydz

Equating this time rate of change of stored energy to the

net rate of energy gain:

3T 3T 3T ‘23 _

‘—— —— + c ——'+ w ) —

0sea at + plcl at O9 9 (3t 32 [9]
3 3T 3 31 _3

The heat flux per unit area due to convection, qcz, at
any point can be expressed as the heat capacity of the fluid

at that temperature times the velocity, or

qcz = ngT

where C q

the convective heat flux with depth 2 is:

- volumetric heat capacity of the gas. The change of

27

3C
3 _ .2! __9 2!
32 (qcz) - ng 32 + WT 82 + CgT 32

Assuming the relative humidity of the air at a constant

100% throughout the depth considered, then 333 g 0

32
Substituting the above expression and corresponding volu-
metric heat capacities (C - PC) for the solid and liquid

constituents brings equation (9) to:

3T 3T 3 3T 3 3T 3T aw
—+ —=—— —— _ — —— _
(Cs+cl+cg) at Cgw z 8x(Axax) + az(xzaz) + ngaz I CgTaz

Letting C - volumetric heat capacity of the composite ele-
ment, substituting the macroscopic gas velocity V for w, and
further assuming an isotropic soil medium (Ax - Az -*), the
heat equation finally becomes:

22 =._1 .12 .2 .12 .2! [101'
C at ax(A ax) + 32‘A 32) + CgT 32

2.5.4 Modification of the Moisture Equation
The infinitesimal element of soil in Figure 2.2 is now
analyzed for mass continuity. The mass flow rate is of the

form:

m = mass of water = 9 water . volume water
areaifime area-time

A volumetric liquid flow rate per area can be expressed
as a diffusivity times a gradient: gradients occur in soil
due to both temperature and moisture. Then the expression

for ix becomes:

3T 39
”2 [DTx 3; I Dex 3x

and the change in mass flow through the element in the x

direction would be:

28

 

 

 

612.6sz azdz , .
T my+amylay-dy
I O
dx
mx >' / {3‘ —-—>rhx+6rinxlax-dx
r'hy

 

 

 

l

rI‘z

Figure 2.2. Infinitesimal Element of Soil for Mass Flow

29

39x 3 3T 39
§§_ dx (dydz) - O2 3; [DTK 5; + Dex 3; dx (dydz)

However, the introduction of air-vapor mass flow
through the soil element due to air density and pressure
gradients produces an additional term in the z direction for

mz:
8T 39
or [DTz 3? + Dez 3'2] + 09. (evv)
and the change in mass flow through the element in the z

direction becomes:

am
2 - a. .32 9.9
—-3—z- dz (dxdy) - pp“ {82 [DTz 32 4» D62 32 + evV]} dz (dxdy)
The air mass flow term is:
‘ a _ aev av
fi(evV)-Vaz +£3sz-

If the moisture content is considered above the wilting
point, the relative humidity of the air-filled pores remains

constant at 100%, and

Assuming an isotropic soil medium where diffusivities are
constant in all directions, the total change in mass flow

through the element is:

a 36 8T 3 3 3T 8V
___. _ + _. + __ __ __ _—
2 3x (D D 3x) 32 (D9 32 + DT az) + 6v 82] dxdydz

p 63x T

The time rate of change of stored mass in the element

is:
___§_ mass of water = __8_ 0 (volume of water)
at in soil element at 2 in element

___ 5%[31(volmne of water )] dxdydz = 0 1g dxdydz

 

volume of soil element 2 at

30

Finally, equating the time rate of change with the

change in mass flow, we obtain the moisture equation:

38 3T

3
6 3x + DT 8x

8

CD

3 36 BT 8V
) + 3? (D6 35 + DT 33) + eV 32 [11]

8
- 3x (D

”-

Mass flow adds a term on the right hand side of this

equation very similar to that added to the heat equation.

2.5.5 Evaluation of the Velocity Gradient. Both equations
(10) and (11) call for the gradient of the air mass velocity
with respect to depth z.aV/Bz. This gradient is now evalu-
ated by differentiating equations (6) and (7) and adding.

For the top few centimeters of soil we can assume the
air permeability and soil porosity changes with depth are
relatively small, and K* and e are then constant. Letting A
- K*A*ao/e and 3* 2 (wt - y*z - 2ny/L*), the partial deriva-

tive of equation [6] becomes:

 

 

 

3V * * * *
-SE = Aexp(-l z)(y (rCosa*-Sin3*)+ A (Cosa*+rSin3 )) [12]
For a similar treatment of equation [7],
V * 2'
du le 3 ( 298 ATAz+Vd1 le)
A2 A2 [13]
*, 2'
- ‘Vgge (AT/Az)+(Vd1/Az) -le/Az
The total velocity gradient is then:
3V * e e . e e a . s
3; = Aexp(-A z)(y (rCosa -S1n3 )+1 (Cosa +r81n3 )) [14]

 

+ i (‘Vgge*(AT/Az) + (Val/A2)2 - le/Az)

31

Equations [10] and [11], with equation [14], now form
the governing equations for moisture-heat interactions in
the seed zone of the soil, accounting for wind gusts and

temperature reversals.

2.6 Summary

This chapter developed the two governing equations for
heat and moisture movement in the soil. A separate term for
vertical air mass movement through the soil was developed
and incorporated into the equations. Chapter 3 further
develops these equations into a set of algebraic equations

suitable for modeling on a computer.

3. MODEL FOR NUMERICAL SOLUTION

3;; Background 9f Finite Element Modeling

The governing equations [10] and [11] are nonlinear,
coupled, second—order partial differential equations in
space and time, and require solution by numerical methods.
Zienkiewicz (1971), Segerlind (1976) and others have shown
the application of the finite element method of solution to
many types of engineering problems governed by equations of
this type.

Typically the governing equation and boundary condi-
tions are incorporated into a functional term. This func-
tional possesses the property that any function which makes
the functional a minimum also satisfies the governing dif-
ferential equation and boundary conditions. Variational
calculus is used for the development. This functional is
given physical meaning in the theory of elasticity when it
is the potential energy of the stressed system, but is less
obvious physically in other engineering problems.

A number of investigators have applied the finite ele-
ment method to moisture and/or heat flow in soils and bio-
logical materials. Comparisons with the finite difference
method have been made recently by Gray and Pinder (1974),
Cushman and Hirkham (1978) and Fall, et al. (1980). Fall,
et al. found for a one-dimensional soil water case the
finite element method was simple, flexible, conditionally
stable but not economical for precision over a long time

period. The finite difference method was found to be

32

33

unconditionally stable,complex, fast, efficient, a computer
space economizer, but could not handle unequal grid sizes.

Treatment of time-dependent problems by the finite
element—functional method reduces a partial differential
equation in space and time to a set of ordinary differential
equations with time as the only independent variable. The
finite difference method is then used to reduce this set of
ordinary differential equations in time to a set of alge-
braic equations readily programmed on the computer. This
approach has been used by many investigators, and most re-
cently by Haghighi and Segerlind (1978), Misra and Young
(1979), Gustafson, et al. (1979), Haghighi (1979), Misra and
Young (1980) and Misra, et al. (1981).

An alternate method has been used lately for solving
finite element problems. The Galerkin weighted residual
method has been used successfully as a more direct way than
the functional method of formulating the governing numerical
equations. This method has the advantage of being able to
treat the time variable in time-dependent problems along
with the dimension variables in the spatial domain within
the same formulation, rather than two separate ones. Thus,
the reduction to a set of algebraic equations for computer
programming is obtained with one formulation, rather than
with two in sequence with each other. This has been
demonstrated by Gray and Pinder (1974), Judah et al. (1975),
Cushman and Kirkham (1978), and Cushman (1979). Other ad-

vantages of the Galerkin method are it retains the solution

34

in the form of functions, retains closer contact with the
physical problem, and provides greater flexibility in the
trial functions.

Solving two or more coupled phenomena simultaneously
has also been done recently. Haghighi and Segerlind (1978)
and Haghighi (1979) studied moisture and heat diffusion and
their effect on stresses in soybeans. Gustafson et al.
(1979) studied temperature and stress interactions in corn
kernels. Misra and Young (1980) and Misra et al. (1981)
studied moisture, shrinkage and stress interactions in soy—
beans.

It is the objective of this research to develop a model
of the heat-moisture interactions in the soil using the

Galerkin method for a two-dimensional, time-dependent case.

212 Mgggl Development

Since the air mass velocity has been determined as an
explicit function of depth and time, equation [14.] can be
evaluated directly at every point in the soil for every time
step using the temperature gradient at that point from the
previous time step. Thus, the velocity gradient BV/az in
the equations [10] and [11] can be evaluated independent of
the moisture-heat interaction. With this approximation, the
problem becomes that of evaluating two field variables,
moisture content and temperature, simultaneously in two
dimensions, x and 2, over the time variable t.

Gray and Pinder (1974) demonstrated the Galerkin method

35

for a two-dimensional, time-dependent case. They used
three-dimensional elements with a simple triangle

representing the x and z coordinates in the spatial domain

 

 

i,j,ksnodes in space domain

0,! =original,final nodes in
time domain

Figure 3.1. Gray and Pinder Model

36

and a third dimension representing the time domain. This
element is shown in Figure 3.1 as a prism of triangular
cross Section.

This element can be thought of as being composed of
three spatial nodes with each node located at two different
positions in time. The resulting three-dimensional space
can be considered as sheets of three-dimensional elements
that lie parallel to the x-z plane and extend over the x-z
domain of the problem. Each sheet has a finite thickness
determined by the time span covered by that sheet. The
problem is then solved sequentially sheet by sheet. Initial
time conditions for one sheet are determined by the solution
generated on the previously considered sheet.

The Galerkin method can be applied to time-dependent
problems in finite element analysis if the unknown field
variable can be approximated by:

1111

$(x,z,t) = Z G.(x,z,t) 4. [15]
i=1 1 1

wherei'E = approximation to unknown field variable in'the
space and time domain
Gi = basis function at node i in the space domain,
dependent on space and time
4 = predetermined value of field variable at node i

i
nn = total nodes in the space-time domain

The basis function in equation [15] is composed of a

spatial and a temporal part. The spatial part is the

37

well-known shape function for a triangular element in a
finite element grid, designated Ni at node i of that
element. The temporal part is designated Bi. Thus the

basis function becomes:

Gi(X.2.t) = Ni(X.2)-Bi(t) [16]

The temporal functional can be expressed by a
Lagrangian polynomial in one dimension, which at node i in

the element has the form:
8 =1? ___—(tr-t)
1 r=l (tr-ti)
rfi
where m - number of time locations for node i. For a linear

time dimension with two time nodes for each spatial node in
the element in Figure 3.1, m = 2, and the temporal function
becomes:

t2_t tl‘t
@ node 11, E—::; @ node 12 [17]

 

B. = _
1 t t1 1

2

In Figure 3.1 the total number of nodes of the element
geometrically is six. However, initial conditions are used
in the time domain in such a way that the known field vari-
able at each node at the end of the previous time period is
assigned that same node at the beginning of the next time
period. Thus, the actual number of unknown nodes to be
determined at each time step for each element equals the
number of spatial nodes per element times (m-l) time nodes.

For the linear time dimension, (m . 2), the calculated nodes

equal three.

38

Galerkin's method obtains an approximate solution to a
differential equation by requiring that the error between
the approximate solution and the true solution be orthogonal
to the functions used in the approximation. Generally, if
we start with a differential equation L(¢) = 0 where L is a
differential operator, and approximate the solution by. I in
equation [15], then the solution L (I) =8 wheree is a resi—
dual or error due to the approximation. In order to mini-
mize g mathematically, we require:

1’ ci e det = o

tR
for each of the basis functions Gi. This integral requires
that 61 must be orthogonal to 8 over the region R over the

given time period t. Thus

If G.L($)det = 0 or "’ GiL($) dxdzdt = 0 [18]
tR 1 tzx

The expressions for L (I) for both moisture and temper-
ature for a single element are obtained by substituting the

following expressions from equation [15]

5(x,z,t) =

i Gi(x,z,t) 91

"mm

1

"MON

T(x,z,t) =

i Gi(x,z,t) Ti

1

into equations [10] and [11].

6
" .. .3 .9. _3_ __3. .32. _3. __3_ .9.
“9) ‘ [8x(DT3x) I 32(DTaz)] i GiTi + [ax(Deax) + azme 32)]

6
3 3V
)3 G.6 - -- (I: G.e ) + e ——
i 1 1 at 1 1 v z
6 6
- a a a 3 3V
= - — _ _ e s + — a e
L(T) [3x()‘3x) + 32(Aaz)]§ GlTl C:g 32 i GlTl

(19)

The moisture expression above assumes the vapor content

6v remains essentially constant over space and time.

A notation more convenient than summation terms is that
of matrices. Thus, we equate for each element:

6

i=1 G191 5 [bioGiijoijGkonf]

fi

CDCDGDCDCDCD

if

jet}
jf

ko
kf)

 

 

l

P-Mm

_ ‘ _ e e
=1 GiTi ‘ [éioGiijoijGkonf:] Tio ‘[G ] {I } [20]
T.
if

 

where the subscripts i, j, and k now represent the three

nodes of the triangle in the space domain, subscripts o and

40

f represent the original and final values of the triangle
moving through the time domain, and the superscript e
denotes the matrix value for a single element in the grid of
the space domain.

Substituting these expressions and thOse of [19] into
equation [18] yields the Galerkin approximation to the soil
heat-moisture interaction for a single element:

eT 8 a a a e e
.1; {law -—<-.—.-—>] M

a 3 a a e e a e e
+ [RUDD 3;) + E—z-(DO-B-E)] [G ] {9 } " 'a-E([G 1 {9; )

+ 9 LV} dxdzdt = o

v 2

III e T 3 3 3 3 e { e}
and tzx [G ] { _3x()‘ —3x) + __an .32) [G I T

+ C §_v_ [Ge] {Te} _ C _3( [Ge] {Te} )} dxdzdt = 0 [21]

The second derivatives in the above equations impose
unnecessary continuity restrictions between elements. These
may be eliminated by the applications of Gauss' theorem.
One of the special forms of Gauss' theorem listed in Potter

and Foss (1975) is:

III $.w dV = IIw-n dA [22]

41

Where W is a vector, A is the surface completely surrounding
the volume V, and h is an outward-pointing unit vector nor-
mal to the elemental area dA. If we let the vector take on

two dimensions with magnitudes as follows:

0)
’U

+ +
1 k

W = M + Q

0202
xlru
0)

Z

where i and K are unit vectors in the x and 2 directions,

then equation [22] can be expanded and rearranged to obtain:

2 2
V 8x 82 V
[23]
+ II (M 22 I 1 Q 32 E). .n dA
A 8x 82

The left hand side of equation [23] now is of the same
form as the three terms in equations [21] which contain the
second derivatives. Using equation [23], these can now be
reduced to first derivatives. For sufficiently small
elements, we can also assume the' diffusivities and
conductivities are constant over each element.

Using the space-time domain De for each element rather
than volume V and Ae to denote the elemental boundary-time
domain, we can evaluate each of the three terms. For

example, the first term becomes:
2

82
III e T 3 e e e T- e
De{[G] UT 3:?- ([G] {T}) + [G] DT a:—-2-([Ge] {T} )}dxdzdt

_ _ III _3, e T 8 e 8 e T 8 e { e}
- 0° DT [3x [G 1 -§; [c l + 3; [G ]='§; [G l T dxdzdt

1’ e T 8 e + + 8 e + + { e}
+ e [G ] DT [%§ [G ]1~n + 3; [G ]k-n T

A

42

Similar expressions can be obtained for the second and third
terms.

If we further let:

.3

at [G]

_ .3 _ _3 _
[GX] - axle]. [Gz] - 32 [G], [Gt] -

then incorporating the above expressions into the equations
[21] and summing over all nt elements in the domain yields

the following:

 

 

Moisture:
2t III 0 ([GelT [Ge]+[Ge]T[Ge]) {we} + o ([GethGe]
e=l De T x x z z 6 x x
e T e e e T e e 3V e T
+ [oz] [62]) {a } + [G 1 [Gt] {9 } - 9V 3; [G ] :]dxdzdt
[24]
- M II ([GelT [Ge] ID {Te} +0 {Ge} I) “’3
eil Ae x T 6 l n
e T e e e + +
+ ([G ] [Gz] (0T {T } + De {9 I )) k-n:]dA = 0
Temperature:
gt III 1([GelT [Ge] + [GelT [Ge]) + CIGelTIGe]
e=l De X X Z Z 1:
3V e T e e
- c9 3; [G 1 [c 1] {T } dxdzdt [25]
nt
- II 6 T e e +.* e T e e +.* =
eil A9 [lIG 1 [ex] {T } 1 n + A [c 1 [oz] {T } k n] dA o

The last summation term in each of the two equations

above represents the boundary values around the x-z domain.

43

The dot products i-E and i°h are the x and z direction
cosines and are unity when the 3 vector is considered along
the x and 2 directions, respectively.

We now return to equations [16] and [20] and delineate

the following:

 

 

 

G 1 _ 3Gio aGif BGjo anf 8Gko aka
[ x - 3x 8x 8x 8x 8x 3x
_ aNi.B aNi.B aNj.B aNj.B aNk.B aNk.B
- 8x 0 8x f 3x 0 8x f 3x 0 3x f
Similarly:

N. N. N. N.
_ 3_£. §_£. 3.1. 3.1. .i_5. 3.5.
[Gz] — [ 32 8o 82 8f 32 8o 32 8f 32 8o 32 8f

If we now identify two terms as found in Segerlind (1976),

N N

ba = 2A 33: and ca = 2A 33% for a = i,j,k and where A = the

area of the triangular element, then:

 

_ —1 O O . . O .
[Gx] ‘ 2A [hi 80 hi Bf bj 8o bj Bf bk 80 bk Sf] [26]

_ _1 I C O O O 0
[G2] ‘ 2A [Ci 80 Ci Bf cj Bo Cj Bf Ck Bo Ck Br] [27]

To delineate the matrix [Gt]' we recognize that the
polynomial governing the movement of the triangular element

through the time domain is of the form:
¢ = 8o¢o + Bf¢f

where 4> denotes the field variable at a particular node

44

passing from one time sheet to the next. From equation

[17], and using local coordinates, the expression becomes:

¢
_ - .2 .E 0
¢ — [(1 At) (At)]{T} [28]
f
where At is the time increment from one sheet to the next.
Then:
dBo _ _._l def = _l .
“HE ' At and ’3? At' and°
[G 1 = —l -N N -N N -N N
t At i i j j k k [29]
Equations [24] and [25] can be further delineated by

carrying out the matrix manipulations within the integrals

using equations [20], [26], [27] and [29].

F'- .-

2 2 2 2 - -
b1 0 bi BoBf bibjBo
T e 1 -
[Ge] [G]=-—- 2 .22--
x x 4A2 bi 08f bi 8f
2 (6X6 matrix)
bjbiBo - -
- [30]
L

 

 

[GelTlGe] -1-—
Z 2 4A
e T e _ _l

[G I [Gt] - At

[GelTlGe] =

 

 

 

45

2 2 2
ci Bo cl BOSE
2
Ci B08f
2
cjcisO
2 2
N1 80 N1 80
--N.2 N.28

(6X6 matrix)

--.. T

-NiNj 80

(6X6 matrix)

 

 

(6X6 matrix)

 

[31]

[32]

[33]

The line integral resulting from the boundary condi-

‘tions term

in each of the equations [24] and [25]

simplified with the following substitutions:

can be

46

%% = [Gx] {a} 3'; + [Gz] {9} fi-K
3% = [GX] {T} 3.3 + [Gz] {y} §.g [34]

where n is the outward normal to the surface.

If we abbreviate the expressions [30] through [33] with

the following:

[GelTlGe] = —i3 [be]. [GelTlGe] - -l—

[CB]
X X 4 z 2 4A

Ge T Ge - —l N , Ge T Ge N
I I [ tI At I 8] I I I I I 88] [35]

and sum over nb boundary elements, the governing equations
finally assume the following integral form:
Moisture:

e=l 2

nt
__L III e e
2 [4A W + I...“ dxdzdt we I. I WIT I)

_lIII {e}_ fl [If eT
+ At tzx [N8] dxdzdt 6 eV 32 tzx [G ] dxdzdt
[36]
nb 11 b b
e T 36 ET ) = 0
bil tL [G 1 det (De 3n + DT an ]

Temperature:

nt
2 ._l§ éé; ([bB] + [c8]) dxdzdt {Te}
e=l 4A

47

_C. III e .. ._ III

+ At tzx [NB] dxdzdt {T } C9 2 tzx [N88]
dxdzdt {Te2] [37]
nb b

- 2 A II Ice]T dIdt ( 33) = o

b=1 tL an

where L . length of the boundary. It is important to note
here that the parameters D , D6 and A in the line integral
representing boundary conditions are not necessarily the
same as those in the surface integral representing condi-

tions within the element.

3;; Integration of Terms in the gods;

3.3.1 Area Integrals. The space domain integrations in
equations [36] and [37] can be carried out using techniques
found in Segerlind (1976). Since the b and c terms in the
first integration are constants and the 8 terms are only

functions of time, then:

1
4A2 éé; ([bB] + [c8] ) dxdzdt = 3% { ([b8] + IcsI) dt

Using the area coordinate system,

 

 

-l "‘ [N8] dxdzdt = A ‘ I -280 280- _80 8° -_-

At tzx 12At t -ng 28f 'Bf ___
-BO 80 -280 —-- dt
b-Bf Bf (6X6 matrix»
J.

 

and

III
tzx

Ice]T dxdzdt = III
tzx

III =
and tzx [N88] dxdzdt

{1’

3.3.2.

"h!

 

[-

2
28° 2808f
2
28f
(symmetrical)

II a
tLlG 1T

Line Integrals.

, II
det tL

 

 

NiBo
NiBf
NjBo
NjBf
Nkso
Nksf

 

L
where ka = Iildj, ij, Lki} .

 

..J

= A I
dxdzdt 3 t
2
BoBf so
2
8f BOBf
2
2808f so
2
28f sosf
2
280
I
det = 2%; t

A

 

 

dt‘

2808f

zefz
.J

 

f .

01 so
01-6f
02-3o
02-8f
03'30

dt

 

@3'Bf

In the line integral above, the boundary element is

assumed small in space and time so that the heat and

moisture flux terms

49

%§ and %% can be assumed constant. The final, generalized

form. above will assume specific values of 0 or 1 for the
.terns Ol, 02 and 03 depending on which of the three sides of

the triangular element is on the boundary.

3.3.3 Time Integrals. Returning to equation [28],

— - .E .E
Bo ' 1 At At
Then, integrating from 0 to At, the length of one time

and 8f =

increment, the integrals become:

II II A.
At 2 4 At 2 At At

I = I at = ‘“ I = A2
080 dt oBf 3 0808f dt 6

 

 

 

 

Thus, the first term in Section 3‘3'1'Z% é ([b8] + [c8]) dt
2 1
= %%K -[BC], where the 2X2 matrix is multiplied by
1 2
every element in the BC matrix, and:
F. 2 1 1P 2 .,
[ bi bibj bibk ci cicj Cick
BC] 2
2 2
b3. b.bk + C3 c ok
2 2
(symmetric) bk (symmetric) ck
t _I L -

The second term becomes:
-1 1 2 l 1

.5 .
24 -1 1 1 2 1

 

 

50

 

 

The third term becomes: .ff
AAt 1
6 41
II
1
III
The fourth term becomes: 2 1 1
2 l l 2 1
A .
—“= [1 2] 1 1 2
72
The boundary term in Section 3.3.2 becomes:
At-Lxx 01]
4 01
02
02['
O3
03

 

I

3.4 Set 2; Algebraic Equations
Substituting the above matrices into equations [36] and,

[37], the governing equations for a single element assume

the following algebraic form:*

*It should be noted here that the individual element equations
are not equalities in the true sense of the word. The conduc-
tance matrix of the temperature term for each element is a con-
tribution to the construction of the global conductance matrix
for the entire grid of elements. The same is true for the
diffusivity matrices and the force vectors. Thus the pseudo
equality symbol * will be used for the element expressions.

 

51

 

 

uoisture:
At 2 1 - [BC] (De{8}+DT{T}) +3-4- '1 1 - 2 1 1{9}
24A 1 2 -1 1 1 2 1
1 1 2
1 01
1 L At 01 *
3V AAt _ xx 18 fl =
" 9V7? 1 "T 02 (De n+DT n) W [38‘
1 02
1 03
1 03
Temperature:
2 1 II
2 1 1 1 1 2 1
lAt . 95 T
""2411 1 2 [BC][T]+24[-1 1 1 1 2 {}
_ _I
2 1 1
2 1
_ g! AAt .
C9 32 72 1 2 1 2 1 [T] [39]
1 1 2
01
AAtL
- xx 222.
4 01 an - {0}
02
02
03
03

 

 

52

These represent a set of six algebraic equations per element

for each of the two field variables, a total system of twelve

coupled equations

nodes, each at two different times.

The above equations can be rearranged to the

form:

Moisture:

D

6

 

 

0)
<1
3’
D
‘1.

<:
2:!
«II

I-

At [
24A [BC]'

 

[£2.23

 

 

 

 

 

 

 

 

 

 

I—II-II—II-Itoro

01
01
02
02
03

03
L

D At
24A

—De

[BC]°

I—II—INMI—II-I

L'"

 

fix

At

 

LNIoI-II-II-II-IJ

0)
CD

02
:3

for moisture and temperature on three

following

io
if
jo
jf

ko

 

 

9
L ka

 

01
01
02
02
03
03

 

0
0.0
mo:
~|<

 

VAA.
7

Mr?

l-‘Nl-‘NNb

 

[—

NHNHDN

HNNnbt-‘N

 

NHONNI-I‘

53

 

ubNNt-‘NH

j

 

 

EH?

io
if
jo
jf
ko
kf

80-30-30-386

 

I

hIId hIIa nauo

AAtL

XX

3n

hIru N «3rd H

8.2

l N «Ira H PJIHJ

 

”II

io
if
jo

Iii-3h]

f

ko
ka

Flt-3P3

 

[0}

These matrices can be combined to produce the governing

equations in the following form:

Moisture:

FT matri
(6 X 6)

U matrix
(6X6)

X

ejO

] {Tic Tif Tjo ij Tko ka}

] I1.

{1 1 1 1 1 1}. T :....._

[40]

AAt
6

as T
o 33){01 01 02 02 03 03}

54

Temperature:

W matrix T
- ['(6 x 6)]{Tio Tif Tjo ij Tko kaI [41]
* At
= 5%— (Igg-I {01 01 02 02 03 03} T

This set of twelve equations has now been given physi-
cal meaning through delineazation in space and time. At
each calculation, six of the field variables are known,
having been solved from the previous time step. Further,
temperature and moisture are now decoupled. Thus, the‘ tem-
perature is solved first, then the moisture can be solved.
This set can be reduced to six equations in six unknowns by.
combining pairs of equations. The first pair of the tem-

perature subset will take on the form:

+WT +WT +WT.+W

W Tio+w T' ”‘7 T' ”'4'” 5 ko 6 kf 7 10 8

1 21f 330 Ti

if f

* ALxxAt 31'
+ W9Tjo+w10Tjf+w11Tko+W12ka = __‘2 SE (01’

where Wi are the elements of the w matrix in equation [41].
The other five pairs of equations assume a similar form.
Rearranging the matrix elements and variables as a
temperature contribution fi‘e) and moisture contribution 134m)
from each element, the governing equations finally assume

the following form:

55

 

 

Temperature:
81 matrix Tif I. At 01 52 matrix T.
(3 X 3) T - _’_‘3.‘.__ (132) + (3 X 3) 10 = R (e)
jf 2 an 02 - T
T f 03 T30
k ko
[42]
$231952:
Rhmitg)” 61f - LxxAt (D 33 + D 33-9) 01
ij 2 T 3n 8 an 02
ekf 03
- e 3V.AAt 1 + R23m2t§1x Bio
v.3; 3 1 ( ) 9'
30
1 6ko
+ [iffy] (:10 : :if: (e) [43]
( j0 if = RM

(
Tko + ka)

The elements of the S temperature matrices each contain
a soil conduction, soil bulk capacitance, and a soil air

convection term characterized by the following example:

5 = AAt(BC11) EA

1 8A 6 g

+
I
0
l2
|>
D
d

0)
N

H
N

The elements of the R moisture matrices include a

liquid diffusivity and a moisture capacitance term of the

56

The elements of the Q moisture matrix contain a thermal

diffusivity term of the form:
Q1 = DT At(::11)

Equations [42] and [43] now comprise the six equations
needed for moisture and temperature at the three nodes of
any given element for a given time.
ggg Summary

This chapter brought the governing equations of Chapter
2 to the algebraic form necessary for modeling on a com-
puter. These equations require the soil parameters A, De ,

D C, Cg, 6v,andaV/az. These parameters are now delineated

T’
in Chapter 4.

4. SOIL PARAMETERS
The soil properties needed for solution of the model
comprised of equations [42] and [43] are D , Dr, A, 6v,
BV/Bz, C, and Cg, all identified earlier. These prOperties,
as they relate generally to all soils or to average temper-

ature and moisture conditions, are now developed in more

detail.

4.1 Moisture Diffusivity QB

Philip and DeVries (1957) showed that De 8 Del + Dev .
However, DeVries (1958) and Jackson, et al. (1975) showed
that DGV' is significant with respect to Del only at
extremely dry soil conditions (see Table 7.2). Thus, De can
be set equal to Del 2 K (aw/391 ) with no loss of generality

in the moisture range above the wilting point.

4.1.1 Moisture Characteristic Slope aw/ael. Ghosh (1980)
proposed a method of calculating soil moisture characteris-
tics from mechanical properties. The method is applicable
for those soils for which the w -e relation can be

expressed as

[44]

where I“; a air entry value of matric potential, mHZO
g) a saturation moisture content = S
8 = dimensionless empirical constant
Ghosh fitted actual data from eight soils ranging from

sand rto clay, and the empirical versus published I values

57

58

gave correlation coefficients ranging from 0.88 to 0.99.
The following relations or values were determined.
6e (air entry water content)1: 0.98, according to Ghosh.

S = total porosity = 1 - pB/(2.65x103)

and

0.2822 0.0625} 0.125

8 = 2.619(12/11) 4

(14+0.7)

0.0625
(5.9113/(11+A3) + 1.1)

where Al 2 3 = % sand, silt, and clay, respectively
I I

A = 6.2 Az/Al - 5.9113/(11 + 13)
The value of we was published for the soils considered.
For any unknown soil, a single w-e reading is sufficient to
establish the relation [44]. Once the relation is estab-
lished, the slope can be obtained mathematically.
~84 58 -84 (1-0 /2650)8
e e B

it _ ___... _

ae e3+1 68+1 [45]

 

4.1.2 Hydraulic Conductivity K. Recent work by Mualem
(1978) and Hirschi and Moore (1980) has established a mathe-
matical relation to obtain K from easily measured soil pro-

perties. The basic relation from Mualem is:

59
where: K - K/K
r 5
KS - K at saturation
Se " 8/80
n - .015w + 3.0 from Mualem (1978)
w = energy required to drain a unit bulk volume of soil
from saturation condition

6
I 0 wad9
6H

w - matric potential, cm 320
Y - weight density of water, g/cm3
9H a moisture content where soil is considered
dry, or where w-e curve reaches the assymptotic
level
Using the expression [44] and integrating:

we 8 -8+1

6 - -
d9 - -s+1‘s S 9H I I451

w = I Y w (--0
H O

‘8

 

and:

4.1.3 hysteresis Effects. Hillel (1977) stated the follow-
ing:

“After its strong daytime desiccation, the surface zone of
the soil draws moisture from below during the night, so that
the top layer is in a process of sorption while the underly-
ing donor layer is in a process of desorption. In princi-
ple, the hysteresis phenomenon makes it possible for a
sorbing zone of soil to approach potential equilibrium with
a desorbing zone of the same soil while the former is at a
lower moisture content, and hence at a lower value of
hydraulic conductivity. It would seem to follow that the
hysteresis effect can contribute to the self-arresting ten-
dency of the evaporation process by causing it to fall below
the potential rate earlier than it would if hysteresis were
nonexistent.”

60

Hillel, in a computer simulation, demonstrated that
hysteresis has a marked effect on the overall evaporation in
a diurnal cycle, and thus cannot be ignored. Bresler et al.
(1969) and others have determined that K is essentially a
unique function of 6, and the hysteretic effects are present
in the w-6 relationship. Therefore, a separate relation,
similar to equation [44] but with different parameters,
needs to be defined and stored for the sorption curve of the
soil in question. Since sorption curves have not received
much attention, are difficult to measure, and transition is
still not well understood, transition between the two curves
when used in the computer simulation will be delayed by one
time step. Further, only the two primary curves will be.

used.

4.2 Moisture Diffusivity g5
Philip and DeVries (1957) showed that QT - ”Tl + DTv'
These diffusivities are further delineated to:

DTl = KYw

DTV = fDatmvBh[VT)a/p1VT

61

The following parameters were either evaluated by
Philip and DeVries or cited by them from earlier works:

h . 1.00 for 4 above the wilting point

0 a 1.024 for T - 20°C

- 1.73 x 10.2 kg m-3 for T - 20°C
3 3 1

po

8 =- 1.05 x 10' kg m“ II c" for 10°C .< T .< 30°C
4.2.1 Flow Factor f._ This factor was equated to:
f - S - 60 for 61 S 61k

f . a + 361/(S-le) . (s-el)(1+61/(s—elk)) for 91 > 91x
since a - S - 61.

4.2.2 Moisture Content elk’ This has been defined as the
point where liquid continuity fails. It is somewhat arbi-
trary, but may be extracted from a K—e performance curve
where K falls to some percentage (perhaps l/1000) of its
saturated value.
Philip and DeVries (1957) cited from Moore's work a elk

- 0.20 for Yolo light clay. DeVries (1958) used 91k . 0.10
for a medium sand. If the point where K - KS/1000 is used,
the corresponding 6 from several actual soil K-e curves
published in Mualem (1978) are tabularized below. The

values in parentheses correspond to e . 40% saturation.

62

Table 4.1. Moisture contents at K a Ks/1000 for
several soils (Mualem (1978))

 

 

Grenville silt loam .34 (.20) Volcanic sand .11
Adelanto loam .22 (.17) Plainfield sand .10
Loamy coarse sand .20 (.16) River sand .08 (.16)
Fine sand .12 (.14) Sand fraction .07 (.14)
Silt Mont Cenis .20 (.18)

 

Table 4.1 demonstrates clearly the wide variability of
moisture parameters in natural soils. Because of this and
the arbitrary nature of liquid continuity, a 91k = 40% of

00, as suggested by Jury and Letey (1979), will be adopted.

4.2.3 Atmospheric Diffusion D Using P = 762 mm Hg

atm'
(standard total pressure for air) and converting T to °C,

the term becomes:

_ -7 2.3 2 -1
Datm,' 5.80 x10 (T + 273) cm sec

01'

2 1

D = 5.80x1511(T + 273)2'3 m sec_ [48]

atm

4.2.4 Air Pore Temperature Gradient Ratio (v T)a /(v T).
Philip and DeVries listed temperature gradient ratios for
different 9 and S values. A linear regression analysis of
their results calculated by the author produced the follow-

ing expression:
(VT)a
(VT)

 

= 0.248 - 1.535 + 2.46 [49]

 

63
This is good for 10°C 5 T s 30°C and 0.1 s e s 0.5.
Their theory produced values for temperature gradient ratios
between 1 and 2.07. However, it should be noted that Hadas
(1968) reported values up to 3.2, Cary (1965) reported up to

5, and Woodside and Kuzmak (1958) reported values up to 20.

4.2.5 Theory Modification. Although some investigators
have tested the Philip-DeVries theory under field conditions
and have proposed correction factors to make the purely
diffusion mechanism better fit convection phenomena, Jury
and Letey (1979) actually studied the theory itself. They
compared the mechanistic theory with that of the thermodyna-
mic approach by Cary and Taylor (1962a, b) to explain vapor
diffusion through porous media under the influence of tem-
perature gradients. These two theories were compared with
experimental data from five investigations.

They explained that discrepancies between the Philip-
DeVries theory and experimental data arise from the "liquid
island” effect, accounted for by the factor f. Whereas
Philip and DeVries considered liquid islands as equivalent
to vapor flow paths and assumed tortuosity was included in
the (IIT)a/(VT) calculation, Jury and Letey argued that
liquid islands are actually shorter paths as seen by diffus-
ing vapor, and the tortuosity needs to be adjusted. They

proposed the following path correction factor:

64

 

a + 619(a) 2
E = a + 61 (1X) [50]
37a) A1
where 9(a) = 1 for e < 91k
= a/(S-elk) for 6 ) elk
Av/Al = 0.124

If the above factor is multiplied with the other fac-
tors in the expression DTV = f Datm v8h (VT)a/01VT, the
theory more closely predicts vapor flow under temperature

gradients.

ngg Thermal Conductivity A;

The thermal conductivity of the soil has been presented
in the literature as a weighted average of the several soil
constituents. A fundamental study spanning many years has
been published by DeVries (1963). The soil conductivity has

been equated to:

_ zxikili

inki

where xi, ki, and Xi are volume fraction, ratio of constitu-

ent temperature gradient to overall gradient, and thermal

conductivity, respectively, of the ith constituent.

65

When a < a then air in the soil pores is considered

lk'
the continuous medium. For a a 61k' the water is considered
the continuous medium. Whichever medium is used as the
basis for 1, the k for that constituent is taken as unity.

Thus:

x A + k x A + k x A
s s s

 

 

av av 1 1 1
A = for e < e [51]
xav + k1x1+ ksxS 1k
and
x111 + kavxavxav + ksxsks
A = x + k x + k x for 6 ; elk [52]
1 av av s s

where subscripts s, l, and av designate soil, liquid water

and air-vapor mixture in the soil pores.

4.3.1 Soil Temperature Gradient Ratios ks. The expression

for ki is given in DeVries (1963) as:

A. -1
k. = 1 Z [1 + (.1 - 1)gN] [s3]

1 3 N=a,b,c 10
where a,b,c are the dimensions of the three axes of the
constituent particle and A0 the thermal conductivity of the
continuous medium.' The quantity ghldepends on the ratio of
the three axes. For soil particles, the dimensions can be
assumed closest to an ellipsoid of revolution where a = b =
nc. The constant of proportionality n can be assumed to be

around 5. From a graphical relation between ga and n shown

66
by DeVries for ellipsoids, ga = 9b = 0.125. And, since qi +
9b + 9c g 1' go a 0075.
The thermal conductivities of several natural elements
are tabularized below from DeVries (1963) and Hadas (1969a,

1977a).

Table 4.2 Thermal conductivities of several elements,

 

 

J/m sec°C
quartz sand 8.799 dry air at 20°C .02577
silt 2.179 water at 20°C .595
clay 2.975 saturated air at 20°C .0997

 

 

 

Using the above values in equation [53], and using

saturated air as the continuous medium for e < .6 the

1k’
resulting ks values are tabularized below.

Table 4.3 Soil temperature gradient ratios

 

 

Element 6 ; elk e < 61k
quartz sand .274 .061
clay .527 .159
silt .611 .205

 

 

 

 

 

4.3.2 Temperature Gradient Ratio of the Dispersed Fluid kav

or k1.

When 6 a k' water is the continuous medium, air is

61
considered a dispersed particle. DeVries showed that as 9

approaches saturation, the shape of air pockets will become

67

spherical, where 9a = 9b a gC = 1/3. When a < elk' air is
the continuous medium, and the water is the dispersed fluid.
As 9 approaches zero, the ga for water becomes approximately
0.013. The actual shape in each case depends on the water
content as:

ga air voids = 0.035 + 0.298 (6/5) for 6 a 81k

Since 9b = 9a and gc = 1 - ga - gb, then from equation [53],

 

 

_ 2 1
kav ’ 0'333[0.971 - 0.248(0/5) + 0.226 + 0.496 (e/sI]

for e a 61k

[54]
Similarly:

ga water = 0.013 + 0.085(6/61k)

and:

 

 

_ .., 2 + l?.........
k1 ‘ 0'3” [1.287 + 1.878(6/61k) 22-514 - 3-755<9/91RJ

[55]

for e < 61k

4.3.3 Thermal Conductivity of Air-filled Pores xav' Philip

and DeVries (1957) proposed:

The vapor term takes into account the vapor distillation due

to temperature gradients in the air-filled pores. And

IV = hvBLDatm [56]

68
All terms on the right hand side have been identified except
the latent heat of vaporization, which is L = 2.451 x 1&5
J/kg at 20°C.

Recent investigators have found this expression under
predicts actual heat flow. Although temperature reversals
and wind gusts under actual field conditions may account for
heat and vapor flow beyond that accountable by the
Philip-DeVries theory, an adjustment to the theory has been
proposed by Hadas (1977) and Sepaskhah and Boersma (1979) to
improve predictions of the diffusion mechanism itself. This

adjustment is the incorporation of a "mass enhancement

factor” § in:

A = A + 51 [57]

Sepaskhah and Boersma suggested§ = 1.3 for loam and silty

clay loam soils.

4.4 Heat Capacities g and 99
DeVries (1963) showed the expression for heat capacity
of soils to be:

c = (1 - S)cS + 0c1 + acg [58]

where C, C8, C1, Cg. are the volumetric heat capacity of the
total soil medium, soil particle constituent, liquid and
air-vapor gas mixture constituents, respectively. The con-

stituent values are tabularized below.

69

. . -3
Table 4.4 Heat capac1t1es, Jm °C

 

Quartz and clay minerals 2.011x106
Organic matter 2.514x106
Water 4.190x106
Air 1.257x103

 

 

 

According to the properties of saturated water vapor
listed by VanWijk (1963) and DeVries, the density of vapor

can be expressed as

1.015 -3 -3

06 = (0.90m )x10 kgm [59]

for 10°C s T s 30°C. Since we can assume the soil air is
saturated over most of the moisture range, and 01: 1000 kgm:3
then the ratio 00/01 = the ratio of liquid volume to total

volume of the mixture.

Then:
c9 = (0.90x10'6 Tl'015)(4,190x106)
+ (l - 0.90x10-6 Tl'015)(1.257xlo3) [60]
3 3°C—1.

At 20°C, c9 = 1.332x10 Jm_

4.5 Vapor Content 8V;

The volumetric vapor contentev_can be obtained directly from

the expression [59].

av =(0.90x10"6 Tl°015)a [61]

70

for 10°C é T s 30°C.

4.6 Velocity Gradient 3V/3z.

 

4:641 Pressure Changes Component 32p; 35; Equation [6]
expresses the velocity component due to pressure fluctua-
tions proposed by Farrell et al. (1966). Both the works of
Farrell and Fukuda (1955) were supported by experimental
data, but both admitted there was insufficient data to go
beyond a qualitative conclusion. For this reason we can
assume a soil of infinite depth with the wavelength of the
air pressure fluctuations of the same magnitude as the soil
thickness. These were assumptions made by Fukuda which make
the analysis simpler than the model by Farrell et a1.
According to Fukuda for equation [6]:

* * * * ’
L =oo,1 =y = We/zxpo

2 t (t = eriod, sec.
n/ p p p )

 

where w
e = a [62]
p0 = 10.35 m H20
* .
Air permeability K is highly variable with porosity,
soil particle size, and moisture content, as shown in the

table below.

71

 

 

Table 4.5 Air permeabilities, K*, m sec-l
Material Particle Air content Permeability Reference]
diameter (mm) a
32:: “"6" 8335?: :32 8:34-33 “15““
:2: 3:322: ::
Sgheres ‘ 16:3 :43 13:33 Fafrell

 

 

 

 

 

 

 

For unstructured soil, the air permeability is much
more sensitive to the smaller particle sizes. With appreci-
able clay and silt, K* can be expected to be quite small.
However, for drier, tilled surface layers with the desired
degree of aggregation of l to 5 mm, K* can vary significant-.
1y.

Pall and Mohsenin (1980) showed data to match a rela-
tionship between permeability and functions of particle
diameter and air porosity. Using this relationship and the
permeability values tabularized above, the following equa-

tion is developed:

* 2 (1.111x106)a3

K = D
vs (1-a)2

[63]

where DV8 - volume surface mean diameter (m), delineated in

Table 4.6 for sand, silt and clay.

72

Table 4.6 Volume surface mean diameters, mm

 

Constituent Particle Size Range Representative Size

 

 

 

 

 

Silt .05'.002 0026
Clay <.002 .001

 

Thus, for an unstructured soil:

3

 

3 3 3 -
= [31.025) 11 + (.026) 12 + (.001) A31x10

D [64]

VS (1.025I7'1l + (.026)212 + (.001)213

Allmaras et a1. (1977), in their study of tillage
effects on soil structure, found tilled layers with a total
porosity greater than 60% (bulk density < 1033 kg/m 3)
exhibited significant air mass turbulent effects in the soil
pores. Ojeniyi and Dexter (1979) studied the size and dis-
tribution of soil aggregates and pores produced by various
tillage treatments. They stated that much of the drying of
tilled soil in the field occurs as a result of convective
transport of air through pores larger than 8 mm. They found
that, for a sandy loam after one pass of a set of tines, the
portion of the total porosity greater than 8 mm in the 10-20
cm. depth may range from 39 to 48%.

Since tillage greatly affects the porosity and air
permeability, K* must be corrected for this condition.
Ojeniyi and Dexter proposed an empirical expression for the

aggregation:

73

B = 25.2 e‘o'192N

where: D a mean aggregate size, mm
N - number of implement passes
For lack of other quantitative expressions, this will be
incorporated into the model for structured soil as:
-0.192N [65]

DVS = 0.0252 e

where: DVS = mean aggregate size, meters, for the tilled

layers.

Making use of the simplifying expressions [62]:

**2
3V 2K A a * *
KER = -——§—-° exp(-A z) Cos(mt-A z) [66]

Farrell et al. made several field measurements of pres-
sure fluctuations at the surface of a grassed field. Probes
were placed at spacings of .05, .15 and .60 m on the surface
and differences in pressure between them recorded for wind
speeds of 4.20 to 6.60 m/sec, measured at a 2-meter height.
Table 4.7 shows several pressure wave values produced by

that study.

74

Table 4.7 Air pressure wave measurements
(Farrell, et al.)

 

 

Wagelength Period Amplitude Frequency
L (m) 't (sec) a (m) w (sec‘ )
p o
m 10.0 .00100 0.628
.60 2.5 .00025 2.513
.30 1.5 .00015 4.189
.15 1.0 .00010 6.283
.05 0.5 .00005 12.566

 

 

 

 

 

 

The term avp/az in the moisture and temperature models
is evaluated at different depths z and time increments. The
time increments chosen are to conveniently characterize the
germination period and in no way attempt to lock into the
sinusoidal times of surface pressure waves. The velocity
gradient is therefore evaluated at the same point in the
pressure wave sinusoidal train at every time increment
throughout the germination period. If this point is chosen
to be where the sinusoidal component is maximum at 1.0, then

equation [66] reduces to:

 

av
__E = - “a .
32 10.35 exP( 20 70K; 2) [67]

4.6.2 Density Changes Component an/az. According to
Hollman (1976), 8* in equation [13] for ideal gases in the

range between 10°C and 30°C, can be equated to:

= .0034 20'1

 

75

Equation [13] becomes:

 

av .
339 = i [VL0667(AT/Az) + (Val/AZ)2 - le/Aa] [68]

Finally equations [67] and [68] are combined:

av _ “a0 ma
3? ’ 10:35 exP“

-———-—;-z) I [69]
20.70K

 

I

. 2
+-1.[ .0667(AT/Az) + (le/Az) - Val/A2]

4.7 Summary

 

Equations [44] through [69] characterize the general
soil parameters needed for the model. These parameters, in
turn, have been developed to now require the following
specific parameters as inputs:

A

A A % sand, silt, and clay

1’ 2' 3

03 bulk density

we air-entry soil moisture matric potential during
desorption

wea air-entry soil moisture matric potential during
adsorption

Ks . hydraulic conductivity at saturation

8H immobile moisture content

w frequency of surface air pressure waves

ao amplitude of surface air pressure waves

In Chapter 5 the meteorological model and the soil

surface interface will be developed.

5. SURFACE CONDITIONS

The soil heat-moisture flow model equations [42] and
[43] include a time-dependent driving force for 6 and T at
the surface throughout the period of study. Providing
realistic surface conditions rather than imposing artificial
conditions will allow the model to be verified with actual
data, thus lending it considerable strength and applica-
bility.

The governing equations can be revised to a form better
suited for analysis at the soil surface. The terms contain-
ing the boundary or surface conditions are:

a. A%% = heat flow to or from the surface of the soil

b. (D 33 + D 32)= evaporation rate E from the surface,

T 6n 6 6n considered negative when the flux is
from the surface into the air. It is
the volume of water per surface area
per time.

eel Evaporation ee £22 Surface

Equation [43] calls for the evaporation rate E at the
surface. Van Bavel and Hillel (1976) and Hillel (1977) have
developed computer models incorporating a straightforward

calculation of the evaporation:

 

s a
m Rc
. -2 -1
where: Em,= mass evaporation rate, kg m sec , (If we
3

divide Em by the density pl = 1000 kg m‘ , the result

is E, m sec-l, the form called for in equation [43].)

76

77

Hs = absolute humidity of air at the soil surface,
kg m-3

Ha = absolute humidity of air, kg m-3

RC = aerodynamic coefficient between surface and

reference elevation

According to Hillel, the absolute humidity of air at

the surface is obtained from the expression:

HS = Ho exp [ms/46.97(TS + 273.16)] [71]

where Ho = saturation absolute humidity at the surface

temperature Ts,

 

and ws = soil matric potential at the surface, m H20. Also:
1.323 exp (17.27T /(237.3 + T ))
H = s S [72]
0 273.16 + T5

The air humidity Ha can be obtained from the relative
humidity h and the air temperature Ta, or from a direct
measurement of the vapor pressure. If we define the rela-
tive humidity h - pv/ps, where pv.= vapor pressure and p3 a
saturated vapor pressure, both in mbar and at temperature Ta,
thenjx-

pv = hps [73]

Using regression analysis, the author formulated an
expression forjpS from a table of values in Mark's Mechani—
cal Engineer's Handbook.

(1.8Ta + 32)

pS = 2.039 (1.0365) mbar [74]

78
The value of pV can then be calculated and substituted

into the following expression for absolute humidity

 

 

0.2175pv
Ha = (Ta+273.16) [75]
The aerodynamic resistance Rc is obtained from:
RC = Ra(St) [76]
where Ra = adiabatic or neutral value of RC, sec 111-1
St = stability correction
The neutral value is given by:
[ 2
1 2.
R = ne( 0/zoiL_ [77]
a 0.168a

where 20 = surface roughness, m
Sa_= wind speed measured at 2m, m sec-l
Both 20 and Sa are inputs to the model. Hillel (1977)
chose a zc>va1ue of .01 m for his model of natural soil.
VanWijk (1963) presented a table of values for a number of
different surface conditions, both bare and grassed, and
Allmaras, et al. (1977) gave values before and after several

tillage treatments, as did Cruse, et al. (1980). Table 5.1

presents several recent 20 studies.

79

Table 5.1 Surface roughness 20, cm.

 

  

Nicollet
clay loam

Condition

Undisturbed
Plowed

Plow -
disk -
harrow

A;

.47

2.92
1.56

     
       

3.47

2.01

Researchers

 

Allmaras, et a1. (1977)

 

 

Input to
odel

 

Bare, undist-

urbed

Plowed

 

2.44

 

Cruse, et a1. (1980)

 

The stability correction St is obtained from:

St

1

=——

(l "' lO-Rl)

where Ri - Richardson Number.

obtained from:

This dimensionless number

9.81(2.0-zo)(Ta-Ts)

 

R1 =

2

(Ta + 273.16)sa

[78]

is

[79]

where Ta - air temperature, °C, an input.

It should be noted in equation [78] that if Ri ;

an absurd answer would result.

This

+ 0.1,

condition would be

indicative of an extreme inversion at very low wind speed.

 

80

To prevent this, the value of R1 was limited in Hillel's

program to +0.08.

5.2 Heat Flux ee EQe Surface

The soil model includes an expression for surface heat
flux, AaT/a n. Heat flows mostly by conduction in the
soil, and by convection and radiation through the air above
the surface. This expression in the model represents sen-
sible heat loss (or gain) by the soil through the surface,
and can be calculated from a surface energy balance analy-
sis. This is stipulated in VanBavel and Hillel (1976) by:
QRN = LEm + QA + QS [80]
where QRN 3 net irradiance, L = latent heat of vaporization,
QA a sensible heat flux to the air, and Q a sensible heat
flux into the soil. S
Equation [80] represents the total heat flow to and from the
surface, and can be used to calculate Qs at each time incre-
ment of the model. The value QS will then be the surface

heat flux appearing as (ABT/Bn) in the model, a positive

value when heat flux is from the surface into the soil.

The net irradiance can be obtained from VanBavel and

Hillel (1976) and Hillel (1977):

_ _ _ 4
QRN — (1 a1)QRG + QRL eo(TS + 273.16) [81]

where QRG - total global short-wave irradiance from the sun
and atmosphere, W m-Z. This can be obtained from meteoro-

logical measurements. Further, QRL.= incoming long-wave

81

radiation from the sky, W m-z, al a albedo or reflectivity
coefficient of the surface, s = emissivity of soil surface,
and 0 - Stefan-Boltzmann constant, 5.67x10”8 W m72°x74.

The albedo of the surface varies with the surface
cover, its color, roughness and moisture content, as shown

in Table 5.2. Van Bavel and Hillel proposed the following

expression for al based on OS of the surface layer:

a1

0.10 for as > 0.25

0.25 for as < 0.10 [82]

0.10 + (0.25 - 68) for 0.10 5 6S g 0.25

Table 5.2. Soil surface albedo

 

 

 

 

 

Condition Albedo Reference
Dark clay, wet .02-.08 VanWijk (1963)
Dark clay, dry .16 "

Sand, wet .09 "

Sand, dry .18 "

Bare fields .12-.25 "

Wet plowed fields .05-.14 "

Heavy clay, 4 nat- .06-.09 Townsend and
ural tillage Migchels (1981)
conditions dur-
ing wheat ger-
mination,

Manitoba

 

 

Based on the Van Wijk values in Table 5.2, ‘expressions
[82] will adequately represent al for 8's considered in the

model.

82
VanWijk further proposed that the earth can be assumed
to be close to a black body in regards to its emission char-
acteristics of long wave radiation. VanBavel and Hillel

proposed the following:

a 0.90 + 0.1808 [83]

where as = moisture content of surface layer
The long-wave radiation from the sky, as developed by

VanBavel and Hillel from a literature search, is:

4
QRL = C(Ta + 273.16) (0.605 + 0.048 1/137o-Ha) [84]

The sensible heat flux to the air is:
Q _ (TS - Ta)ca
A ‘ RC
where Ca. = volumetric heat capacity of air. If we

 

substitute Ca.= 1.257x103 J m'3IIC'1 , then:
QA = 1257 (TS - Ta)/RC [85]

5.3 Summary
Equations [70] through [85] represent the general rela-
tions of a soil-atmosphere model. These equations have been

delineated to require the following specific inputs:

h air relative humidity
Ta air temperature
Sa air wind speed

QRG total global short wave irradiance

2 surface roughness

O

83
Albedo was incorporated into the model as a function of
the moisture content of the soil surface.
Chapter 6 explains the computer model and the

incorporation of equations [42] through [85] into the model.

6. COMPUTER IMPLEMENTATION

6.1 Finite Element Grid of Soil Furrow Profile.

The furrow profile modeled was assumed to characterize
a field of soybeans planted in rows of 70 cm (28 inches)
spacing and at a depth of 5 cm (2 inches). The study of
moisture and heat flow had to be extended below the seed
depth to a depth where changes are no longer significant due
to diurnal surface fluctuations throughout the germination/
emergence period. Based on the investigations tabularized
below and on the model verification study, this depth of
influence was set at 30 cm. This provided a compromise
between increasing programming complexity for larger grids.

and decreasing accuracy with a grid too shallow.

Table 6.1. Depth of influence

 

InvestigatorI Comment on depth of influence

 

Hadas (1968) 6 to 7 cm. from source of change. Con-
sidered a 20 cm soil column semi-infinite.

Rose (1968a) 10 cm for 6 variations, but greater than
15 cm for T variations. States that daily
T fluctuations go to 30 cm in a bare soil.

Hadas (1969a) 6 cm from heat source in lab study.
Kimball, et al. at 5 cm, diurnal amplitude of T and 6 = 73%
(1976) of surface. At 20 am, it was 20%.

Bruce, et al. 13 cm for 6 in tilled soil. 12 cm for T in
(1977) tilled soil. Some water coming from below

15 cm to supply evaporation. Appreciable
change in T and 6 at 14 cm within 2 days.

Hadas (1977a) r20 cm for 6, greater than 18 cm for T.

[Merva (1975) from 5 to 1] cm for daily T fluctuations.

 

 

 

 

' 84

85

The finite element domain and grid are shown in Figures
6.1 through 6.3. The validation and simulation grids were
similar. The simulation grid was designed for more detail
and accuracy in the vicinity of the seed near the upper
right hand corner, whereas the validation grid was used

exclusively for one-dimensional flow.

  

70 cm I” 5cm

 

   

/////<\\\\
————— I -——--—II-----—-——-— 4“]
I : I: |1130
‘ 5 ‘1“
finite element
domain

Figure 6.1. Finite Element Domain of Soil Profile

86

48 elements
36 nodes

“3
38
_s
h)
0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

F-==_., ._:=-1!-=:_.
('3‘ 241

3f) 33:

. 6

1O

20

30

Figure 6.2. Finite Element Grid for Validation

87

48 elements

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

35 nodes
cm.
35 15 4 2 0
_:35a:--- 1
if 2
‘_fl==!:=5-— ‘3
l 5
1O
20
BO

Figure 6.3. Finite Element Grid for Simulation

88

6.2 Boundary and Initial Conditions
Boundary and initial conditions for the finite element

domain are listed in Table 6.2

Table 6.2. Boundary and initial conditions for the model.

 

 

 

Location Condition Application
Each tilled node 9 = 9i, T = Ti for t = 0
in grid
Each untilled node assumed gradient* for t = 0
in grid
Both side and BG/Bn = 0
bottom boundaries aT/an = 0 for all t
Top boundary time-varying heat for all t
flux and evapora-
tion (Chap. 5) ‘

 

 

 

 

*Gradient based on Bruce, et al. (1977) and trial solutions
of the computer model.

There will normally be a gradient of both 8 and T
throughout the soil profile prior to planting. The model
assumes that the planting/tillage operation thoroughly mixes
the soil so that the tilled soil profile immediately after
the planter has passed has a uniform 8 and T. Untilled soil
in a no-till pass retains its original gradient. The values

9.

1 and Ti for all tilled nodes in the grid, including sur-

face nodes, assume a value equal to the average between that
at the surface and bottom boundaries of the grid prior to

planting. This provides a surface temperature as an initial

89
condition used in the heat balance analysis of the surface-

atmosphere interface.

6.3 Specific Inputs 32 the Model.
The summary sections 4.7 and 5.3 list the specific
inputs to the computer model. The time-varying parameters

h, T Sa and QRG are modeled as sinusoidal expressions to

a.
represent a diurnal cycle. Refer to Appendix D for these
expressions. Although the computer model makes rather com-
plex and intricate calculations at each step, the inputs as
listed above are generally easily measured in the field or

lab, and are supported by an abundance of published data.

6.4 Computer :12! Diagram.

The computer flow diagram used to develop the program
is listed in Appendix A. Briefly, the program follows this
sequence:

a. Reads in constant soil properties and administra-
tive limits.

b. Establishes nodes and elements in the grid.

c. For first time step, and using initial conditions
at each node and time-varying conditions at the
surface, calculates temperatures at each node in
the grid, solving total matrix.

d. For first time step, using initial conditions,
surface conditions, and temperatures calculated in
c., calculates moistures at each node in the grid,
solving total matrix.

e. Repeats c. and d.for each successive time step,
using calculated T's and 8's of previous step as
initial conditions of succeeding step, and using
surface conditions programmed for that particular
step.

90

6.5 Computer Program Listing

The computer program and its documentation are listed
in Appendix B. It was written in the FORTRAN IV language
for the IBM 360/40 system at Delta College. Although this
system was readily accessible by the author, it did pose
limitations. Memory capability did not allow a grid finer
than about 40 nodes, and length of runs, because of other
demands of computer usage, were limited mostly to 3-day
simulations.

Six subroutines are used to calculate soil parameters

or surface conditions. These are:

DIFFM produces K, D8 and 6
ll

DIFFT DT

EVAP " E, QS

VEL " 6V/az

THCON " A

' I!

CAP C, 6v

Two other subroutines, DCMPBD and SLVBD, earlier developed
by Segerlind (1976), were used to decompose and solve the
simultaneous equations in the global matrix.

A special feature of the program is the incorporation
of a hysteresis effect. An adsorption moisture curve is
entered as an adsorption we in a DATA statement. Then in
DIFFM a decision is made on which we to use, depending on

whether that element is losing or gaining moisture.

91
ggg Summary.

This chapter explained the development of the computer
program to solve soil water and temperature interactions.
Special programs, separate from the main program, were
developed to study the sensitivity of key parameters.

Results of these studies are given in Chapter 7.

7. PARAMETER SENSITIVITY ANALYSIS
The soil parameters and surface conditions were ana-
lyzed individually using their respective subroutines.
Computer values were compared with published experimental

values.

7.1 Soil Moisture Diffusivity 20L Matric Potential $2.229

Hydraulic Conductivity 5;

Figures 7.1 through 7.3 compare the output from DIFFM
with the experimental values of Adelanto loam from Jackson
(1973). This soil did not exhibit a sharp break in the
moisture curve, thus the air-entry matric potential we had
to be estimated. This estimate was made over a relatively
wide range between -.05 to -.40 m H20.

Several computer runs were made at different We values.
Three examples are plotted. The plots clearly show how
sensitive De, I and K are to the we value, and further indi-
cate an optimum wetexists for the adequate modeling of all
three parameters. That optimum for Adelanto loam appeared
to be around (@3- -.1017 m H20. Higher magnitudes would
improve the 4-6 model but reduce the De model accuracy, and
vice versa.

The relationships depicted in Figures 7.1 through 7.3
are consistent with the state of the art of soil modeling as
reported by Mualem (1978). He proposed the K-8 relationship
which improved predictions over two earlier models by
several orders of magnitude. Yet it deviated from experi-

mental curves of some soils by as much as two orders of

92

 

 

 

 

 

 

 

\ \
\ \
\ \ \
\ \ \
\ \ \\
\
0. \. \ \
9 A
' \ \ 9°. \
O. \ \x 70\X
0“ If.) \ \ o, \
:I: \ \ \\
E \\ \ \
‘ \
3..
To. \\ \\ \\
'2‘ o \ 3%.: \
2 ‘T \ \‘olo
o \ \b
O. K) \\_ 2\\€?
. W \
3 C.) —-model \Qe‘ \ \
‘5' -—Jackson’s loam '03 \\
2 \‘99 \
<2
\
\
8 . . . _r
I.1O .15 .20 .25 .30
Moisture Contente. m’lm’
Figure 7.1. Moisture Characteristics Curve for

Adelanto Loam

94

V

.30 .35

 

 

 

 

dfllv /V [ARON-1“ W‘J 4‘0
I /l. I I! DOE
“U.3
n4
.elc
e 3
Gary
nu.)

[INII «WAGS Omxumfi

 

 

 

 

 

 

 

Moisture Content , e, m’lm’

 

 

602E .v..>:>zo:ocou 2.39.6»...

O
l I]. 2.
Illlllllllliiilll pk. F .1" .nlll.li.lila.llil
[11”] 0.. a? 3.00:: (Lb.
Emmmo.-... 647.. ooE ] rm
Baum Toe Toe . 9.0.
N O_mum 6710—. S [OP «Aloe

Figure 7.2.

Hydraulic Conductivity for

Adelanto Loam

95

 

10"7

, m’lsec

10'“I

 

 

Moisture Diffusivity , D

 

 

 

3"
9 V I V 1 j
.10 .15 .20 .25 .30 .35
Moisture Content, 9, m’lm’
Figure 7.3. Moisture Diffusivity D6 for Adelanto Loam

96

magnitude. The relation in expression [44] itself is not
universal and is an inadequate model for some soils,
according to Ghosh (1980).

It must be further emphasized that D6 is very sensitive

to several properties which are either difficult to measure

accurately or vary considerably with even minute changes in
moisture content.

The hydraulic conductivity changes by perhaps six to
eight orders of magnitude from saturation to wilting point.
It is very sensitive to the exponent in expression [47].
This exponent, in turn, is primarily dependent on w, the
energy required to completely drain the soil. This "dried“
condition and the Corresponding 8H itself is nebulous.
Mualem calls this point the immobile or residual moisture
content. A check with Jackson's Adelanto loam indicates
this points falls midway between the wilting point and the
point where the soil specific surface is covered by a one-
molecule layer of water. The values for I change by two
orders of magnitude within this range.

Since K is strongly dependent on the w—6 curve, and
dee is directly related, De is very sensitive to the assumed
relation between 6 and 6. This relation depends oni3 and
we. The variable 8 is a function only of soil texture, which
can be measured accurately in the lab.

Soil moisture properties are also very sensitive to the

bulk density. For example, according to the model £0t

97
Adelanto loam at 15% moisture, increasing the bulk density
from 1100 to 1600 kg/m3 increased the (an expression in
equation [47] by 5 orders of magnitude.
Table 7.1 displays the sensitivity of the soil moisture
properties using soil #4005 of MUalem (1978) which is

apparently the soil used by Bruce, et al. (1977).

Calculations are based on the expressions used in the

computer model.

Table 7.1. Soil moisture properties.

 

 

we (m) 8H Bulk gensity Exponent K (m/sec) D 6
(kg/m ) n 8 6 =.02
(Mualem) (mzsec-l)
o.1017 .025 1200 8.1 4.63-17 1.9E-12
.020 " 10.0 8.6E-20
-.2000 .323 " 13.1 3.0E-24
" . n " 8.0 6.4E-17 5.3E-12
1100 8.9 1.8E-18 1.8E-13

 

 

 

 

 

 

 

 

At low moisture contents, Table 7.1 demonstrates the
keen sensitivity of K and De with we and 83, both of which

are subjective observations obtained from the v-O curve.

7.2. Soil Moisture Diffusivity QT;
Table 7.2 shows the individual diffusivity components
for two different soils found in the literature. The num-

bers are rounded off for convenient gross comparisons. Both

98

soils demonstrate that DTV is significant compared to DTl]
only at low moistures, generally at or below the wilting
point for those soils (14% for Jackson's loam). Otherwise,
' throughout the useful range of moisture content, its contri-
bution is negligible.

The computer model for DT for Adelanto loam remained
relatively constant throughout the moisture range studied.
Also, comparing Dq.values in Table 7.2 with De values in
Figure 7.3 for the same soil modeled on the computer, the
model demonstrates that D6 is three to.four orders of magni-
tude smaller than 96' This conforms with published state-
ments that the moisture diffusivity contribution due to a
temperature gradient will be small in soils with adequate“
moisture and moderate temperatures.

2

Table 7.2. Moisture diffusivities, m sec-l/(°C or 8)

 

 

 

 

 

 

DeVries "Medium Jackson loam
sand" @ 20’C @ 25°C
D01 DTl Orv Del Dev DTl DTv DT

.05 18-5 ' 38-10 28-11 68-10 28-9 78-11 78-11
.10 " 88-9 ” 18-9 58-10 88-11 88-11
.15 28-5 28-8 " 28-8 18-10 28-12 58-11 58-11
.20 38-5 68-8 18-11 28-7 28-10 18-11 28-10
.25 18-6 18-10 38-12 18-10
.30 68-5 28-7 18-11
.35 168-5 " ”

 

 

 

 

 

 

 

 

99

122; e911 Thermal Conductivity A;

Figures 7.4 and 7.5 display A vs. 8 at two different
temperatures. DeVries (1963) plotted experimental and cal-
culated values for a quartz sand, and these are replotted.
With some refinements of his expressions (equations [73] to
[75]), the model based on water medium came close to his
calculated at 40°C, and in fact improved on his at 20°C.

DeVries gave a convincing argument that at very low
moisture contents saturated air rather than water becomes
the continuous medium for thermal conductivity. At
extremely low moisture dry air can be considered the medium.
Which medium is considered makes a significant difference in
1. Further, at what moisture content this should occur is
not clear and is somewhat arbitrary.

The plots show model studies with saturated air as the
medium at 8 < 81k. These values, based on expressions [51]
and [55], were considerably lower. DeVries also experienced
difficulty with modeling dry soils, and resorted to a multi—
plier of 1.25 to bring theory in line with experimental
values. However, when basing the calculations on water
medium throughout the useful moisture range, the author
produced a better model, which was then adopted for the

final computer simulations.

2,0 2.5

1;.

Thermal Conductivity, A , JImsec-‘C
1.0

0.5

 

100

DeVfiesUQGB)
experimental

DeVries calculated

/ model based
on water medium

,’ \— model at e<elk based
on saturated air

assumed elk

 

c>0 .65 .1'0 .125 * .2'0 .2'5 .3'0
Moisture Content,e, m3] m3
Figure 7.4. Thermal Conductivity of Quartz Sand

8 20°C

101

DeVflesUQSB)
experimental

    
 

 

 

 

 

DeVries calculated

.9
U
4:
II.)
E 33+ -
g
60
£903“ //
.2
13 /’//\-model based on
3 K2,, / water medium
E 1—
8
E 8‘ "7 Kmodel at e<elk
E based on
0 saturated air
5 n.

0'

assumed eI
L k
0 v

o .05 .1'0 .15 ' .20 .2'5 .30
Moisture Content,6, m’/ m3

Figure 7.5. Thermal Conductivity of Quartz Sand
@ 40°C

102

112; §2£l.§22£ Capacity 9;
Modeling DeVries' "medium sand" produced C values from
1.3686 to 2.8336 J/m3°C. for 8 between .05 and .40. This

-compares favorably with the following DeVries values:

Graded Ottawa sand 1.8286 J/m3°c
Crushed quartz 2.1136 J/m3°€
Northway fine sand 2.2836 J/m3°C
Quartz and clay minerals 2-0136 J/m3°C

7.5. Velocity Gradients av /35 and 328135 f Soil Air.

 

 

Velocity gradients due to. pressure gradients were
modeled for a fictitious soil approximating the air porosity
used by Farrell (1966). Comparisons are made in Table 7.3,
showing very close agreement.

Velocity gradients due to positive density gradients
were modeled and proved to be consistently larger than those
due to pressure gradients by several orders of magnitude.
Table 7.3 compares the two at several depths, maintaining a
constant porosity, and using the temperature gradient
reported by Bruce et al. (1977).

Pressure gradients prevail during the day when density
gradients in the soil are negative, but density gradients
prevail at night when the wind has quieted down. Neverthe-
less, moisture vapor movement due to pressure gradients is
expected to be small compared to that due to positive

density gradients.

103

 

 

 

 

 

 

 

 

 

Table 7.3. Velocity gradients, m/sec-m
an/az comparisons at z=5 cm.
Amp. Freq. 8Vp/az avp/az
model Farrell
1.0E-3 .628 5.98E-5 5.94E-5
2.5E-4 2.513 5.89E-5 5.808-5
1.0E-4 6.283 5.79E-5 5.65E-S
Depth comparisons
z(cm) AT avd/az avp/az
0 - l 2 3.85 6.05E-5
l - 2 0 0 6.04E-5
2 - 3 .5 1.93 6.02E-5
3 - 4 .5 1.93 6.018-5
4 - 6 l 1.93 5.988-5

 

 

 

 

 

 

119; Evaperation egg geil Heep Qeig 2; eye Surface.
Computer simulations were made using the published
atmospheric and soil data of VanBavel and Hillel (1976) for
a loam soil in Binghamton, NY. Total global irradiance, air
temperature and relative humidity were entered as sinusoidal
functions of time throughout the day* using, in some cases,
only daily or even monthly averages found in the publication
for experimental verification. Figure 7.6 shows the evapor-
ation model peaked at the same level as the simulation pub-

lished by VanBavel and Hillel, but was delayed by about 2

 

*Refer to Appendix D

104

10'7
\
\
I
I
/
/
/

 

 

 

 

 

 

 

/I \ \
/ \\
/
/
/
/
/
/
/

u /
8 /

Q
5:0 A
c s- I
g / model with air
3 / temp. =16°C am.
0 /
CL
m / l
0’] /

VanBave|& Hillel(1976)
”typical" result 5

‘i‘

O. . , - , .

" M4 8 12 16 20

Diurnal Standard Time,hrs.

Figure 7.6. Daily Evaporation for a Loam in
Binghamton, NY in June

105

hours. In the afternoon, evaporation in the model declined
more slowly. This cannot be discounted as error in the
model, however, since a direct cause-effect comparison was
’ impossible.

It must also be emphasized that certain soil parameters

were held artificially constant during this initial

verification check. These were the soil surface moisture
and matric potential, and the wind speed. The combination of
air and soil surface temperatures used provided no inver-
sion. Had the air temperature been allowed to exceed the
surface temperature, an inversion would have resulted,
making Ri positive. If the wind speed would be decreased,
as at night, Ri would increase more positively, increasing
St and consequently RC, and finally reducing E throughout
the nighttime.

A check of surface roughness and wind speed effects on
the aerodynamic reSistance RC, using Ts—Ta - 4°C, a moderate
differential, produced the values in Table 7.4. For this

combination, wind speed had negligible effect.

Table 7.4 Aerodynamic resistance for AT = 4°C.

 

 

 

 

 

-1 . -l
z°(m) Sa(m sec ) R1 Rc(sec m )
.010 1.22 -.173 52.7
.025 " -.172 36.2
.035 " -.l71 30.9
.025 2.00 -.064 36.6

 

 

 

 

8. MODEL VALIDATION

The computer model was validated using the field data
published by Bruce, et al. (1977). Additional information
on the soil used by Bruce was obtained from Harper, et al.

(1976).

8_.1. Lemma

The Bruce study involved one of the most detailed and
complete works found in the literature on actual soil mois-
ture and temperature measurements. The study spanned eight
days, from June 16 to June 23, measuring moisture hourly at
ten depth intervals to a maximum of 15 cm., and temperature
every 15 minutes at eight depth intervals, also to a maximum
of 15 centimeters.

Soil conditions published were:

1. texture and particle size distribution measured to

2. iglfimdensity vs. depth to 15 cm.

3. moisture characteristic curve

4. saturated hydraulic conductivity

This provided all the soil inputs required by the model
except the adsorption curve for hysteresis. For the valida-
tion, hysteresis was assumed negligible. The author felt
this to be the better option rather than assuming some un-

substantiated relation.

Above-surface conditions published in the Bruce study

were:
1. air vapor pressure
2. air temperature
3. incoming solar radiation
4. wind speed

106

107

Relative humidity for the model had to be calculated
from the vapor pressure using Equation [73]) and surface
roughness had to be estimated. This fulfilled the surface
condition inputs to the model.

A particular day in the Bruce study, June 18, was
selected for validation. This day followed several unevent-
ful days of drying and preceeded a heavy rain by one day. A
copy of the raw data for that day was obtained from Bruce

for better accuracy.

8.2 Mbdel Refinements

The modeled 24-hour histories of temperature and mois-
ture for June 18 are shown in Figures 8.1 and 8.2. The
model accepted the initial conditions shown on the 0000 hr
(midnight) curve, used the finite element grid shown in
Figure 6.2, was exposed to the sinusoidal surface conditions
and stepped through the day in 30-minute increments.

Validation of the model with the Bruce data required
some refinements of the model. These refinements concen-
trated mainly on five areas, although all elements were
eventually scrutinized before being accepted for the model.
These five area were:

1. oscillations of the finite element formulation

2. the air mass flow in the soil

3. the surface conditions submodel

4. filling in where the Bruce data was spotty
5. moisture diffusivity due to vapor

108

Temperature , °C

 

 

 

 

 

 

 

 

 

 

20 24 28 32 . 36 1 4O 44
o. - "w‘ t—
. 9‘6 ‘9 ’
N- . ‘ f
v o , ,
e. . $8 III ..
l\ O , I
l O ,’ I
O ‘ g 4 I
s— o" I f
I l I H
1 \ CNNFS
. l‘ I ,’ ,’ . 0600
E , ‘I I - 1200
0. \ I“ I . 1800
5 i ll I, I‘ 2400
o. l I] \ ,
80 ] I, ‘l [I
N I l
I ' l l
l ,’ I; —— Bruce data
I , l, ---model
I , [
l I l
' l (l
:I I]
g sJ———Aa':

 

 

Figure 8.1.

 

June 18 Temperatures Comparing Model
with Bruce Data

109

Moisture Content ,9, m’/ m’
0 .04 ”018

1 .1? .16 . .29 ‘ .2‘4

      
  

 

 

 

 

 

 

 

 

E
U
E.
4.:
a
0
C)
O
N
—— Bruce data "I
— -— model 3;
III
III
fll
fll
I:
I
8 ~—

Figure 8.2. June 18 Moistures Comparing Model
with Bruce Data

 

110
8.2.1 Finite Element Formulation. The surface conditions
(soil-air interface) model imposed upon the soil model pro-
duced unacceptable swings in the surface temperatures and
moistures. Several different grid and time step combina-
tions were employed, as well as changing the element capaci-
tance matrix to the non-consistent (lumped) form. Oscilla-
tions were virtually eliminated when the time increment was
reduced to 10 minutes, but this would require too much com-
puter time for a 7— to lO-day simulation. A compromise was
reached by using a 30-minute time increment and limiting the
temperature change at any node to 4°C for each time incre-
ment. This produced sufficient stability for the model with
a reasonable physical constraint, based on field data of

soil surface temperatures.

8.2.2 Soil Air Mass Flow. Vertical air mass movement in
the soil was incorporated into the mathematical equations
for both temperature and moisture. This term, including
aV/az, represented moisture movement towards the surface
during the evening and early morning hours due to density
gradients, and a much smaller component due to pressure gra-
dients. The comparable temperature term represented the
corresponding thermal convection due to this air mass move—
ment.

Use of these terms in the model produced moistures and
temperatures higher than those published by Bruce. Compari-

sons with field data indicated that the "mass enhancement

lll

factor” already added to the soil thermal conductivity term
apparently fully accounted for any thermal convection.
Thus, the extra 3V/az term was deleted for temperature. The
density gradient component of the term in the moisture equa-
tion was adjusted to coincide with field data. The pressure

gradient component was left unchanged.

8.2.3 Surface Conditions Model. Because of the lower
inertia of the properties of the air above the soil and the
greater sensitivity of those properties to temperature, this
model experienced greater oscillations than the soil model.
Although these were substantially attenuated by controlling
the surface temperature (see Section 8.2.1), additional
controls were necessary. Limits were placed on the sensible
heat flux to the air and on the evaporation consistent with
maximum and minimum measured values for similar regions

published by VanBavel and Hillel (1976).

8.2.4 Bruce Data. Although the Bruce study was the most
complete available, some information was either missing or
sketchy. Data on the soil moisture characteristics, bulk
density, and soil texture in the sublayer, as well as the
presentation of the temperature and moisture values them-
selves still posed cause-effect uncertainties and required
clarification.

a. §gil moisture characteristics. The published w—e—
curve, measured from -.01 to -15 bar, did not exhibit a

clearly defined IQ! value, thus one had to be estimated.

 

112

Extrapolation of the curve for on was virtually impossible,
thus a computed value for other sandy loam soils ‘published
by Mualem (1978) was used.

The saturated hydraulic conductivity KS supplied by
Bruce seemed to produce simulated moisture in the upper soil
layers higher than published. Comparisons with the texture
triangle and with Ks values published by Rawls, et al.
(1982) showed that the texture and wilting point moisture
content of Bruce's Cecil sandy loam indeed placed it
squarely in the realm of a sandy loam, but the Ks value 'was
far too high, characteristic more of a loamy sand. The K;
characterizing a sandy loam published by Rawls was then used
for the computer model.

Although the data used were that of June 18 in the
Bruce study, the Ire relationship was not measured until
October 30. Bulk densities in the top soil layers had in-
creased during that time period.

b. 22;! density. The bulk density of the Bruce study-
exhibited the high degree of variability in the top several
centimeters of soil and the probable difficulty of sampling
near the surface, especially in tilled soil. At a 1.5 cm
depth, the bulk density, 1 one standard deviation, ranged
from 1060 to 1280 kg/m3 . According to Table 7.1, this
would correspond to a change in D6 of the order of one

magnitude.

c. The g1 horizon. Soil texture of Bruce's Cecil

loam was sampled at several different depths to a total depth

113

of 15 cm The texture was constant down to 9 cm., and be-
came slightly sandier at the 15 cm. depth. Initially, the
texture published for the top nine cm was used for the
-entire 30 cm. profile simulated on the computer. Computed
moistures were typically higher than published values.
However, an earlier study, Bruce (1972), showed a El
horizon at 18 to 33 cm. depth to be a sandy clay loam under-

llying an A horizon of loamy coarse sand, all typified by a

P
Cecil loamy sand profile. Although the 81 texture of the
site in the study may not have been the same as the 1977
study, incorporation of it into the model substantially
improved the moisture values.

d. Moisture and temperature measurements. Soil water'
content was averaged by Bruce over five sampling sites and
temperature over two sites. Thus, a direct cause-effect
relation using one site and its corresponding bulk density,
etc. was not established. Difficulty in sampling moisture
near the surface required extracting samples above 2 cm.
with a spoon. Further, the plotted moisture data had been

Fourier transform smoothed, thus the published plot might

appear smoother than actual conditions.

8.2.5 Moisture Diffusivity Due to Vapor.

In section 4.1 the assumption was made that since
moisture flow due to vapor was insignificant compared to
liquid flow at soil moistures above the wilting point, then

De was set equal to Del' With this diffusivity, the model

114

exhibited moistures in the upper soil layers higher than
actual field data from Bruce. Moistures in the top one
centimeter, according to Bruce, were reduced to below an .
In these top layers, then, liquid movement essentially
stopped, and the model accumulated moisture just below the
surface from moisture movement from lower layers. Recon-
sideration of the basis of the original assumption indicated
that in the top layers vapor movement is apparently a signi-
ficant component when the surface dries, and cannot be
ignored.

A comparison between D) based on liquid movement ionly
and that with vapor movement added is shown in Figure 8.3.
Ignoring vapor flow, the model deviated from the Bruce data
by about 3% at the S-cm. depth, but was within 1 1/2% when
vapor was included.

Thus, the moisture diffusivity DS was reinstated in the

model as

De = Del + Dev
where Dev a a aDatm vgpv (3¢/39)/01RT
With this addition of the vapor component, the model
came into close agreement with the Bruce data as shown in

Figure 8.2.

8;; Comparisons Between the Model and Published Data.

After incorporating the refinements discussed in the
previous section, the model followed actual data reasonably
well. Some probable reasons for the discrepancies between

the model and actual data have already been discussed.

115

Moisture Content,6 , m’lm‘
0 .04 .08 .12 .10 L30 _ .24

  

 

 

 

 

 

 

 

O A
(\l
V
go.
no.
9
. \
E \‘
I; \)
a \‘
I
8 I
O ‘.
(V [I
Bruce :
——---model De‘DeI :
I
I
I
o i
00 V
Figure 8.3.

June 18 Noon Moistures Comparing Liquid
and Vapor Diffusivities

116

In Figure 8.1 the modeled temperature at the surface
did not reach the daytime levels of Bruce, but in the 1- to
6-cm. depth the modeled noon temperature was 2 to 3° C
higher. Generally, the total range throughout the day was
within 3" of the Bruce range throughout most of the soil
profile compared. This model accuracy is consistent with
the objective sought by Cruse, et al. (1980) which was to
develop a soil temperature model with an accuracy of 1 2° C
at the 5-cm depth.

The moisture comparisons in Figure 8.2 indicate. agree-
ment between the model and the Bruce data. The Bruce data
exhibited a greater moisture rise due to density gradients
at the surface between midnight and 6 a.m., whereas the
model exhibited a greater rise from 6 p.m. to midnight.
Model moistures were within 1% of Bruce moistures at the 10-

cm and 6-cm depths, and within 2% at the 3-cm depth.

8.4 Summary
This chapter compared the model with field data

reported by Bruce, et al. (1977). With some refinements,
the model demonstrated good agreement with the data, as
evidence by Figures 8.1 and 8.2. Versatility of the model
will be demonstrated in Chapter 9 through simulations of

different conditions.

9. SIMULATIONS
One of the principle uses of a model is simulation.
Chapter 1 lists several potential applications of a simula-
tion tool such as this. This chapter demonstrates the
utility and versatility of the model for Michigan soils and

conditions.

9.1 Model Parameters
9.1.1 Soil Profile

The soil profile used in the simulations reported here
was previously explained in Chapter 6. The finite element
grid shown in Figure 6.3,isIaslice through the furrow
extending from the furrow centerline out to the midpoint
between rows, and down to a 30-cm. depth. The time incre-
ment for successive calculations of moisture and temperature
was set for 30 minutes, for a total span of three days. All
simulations were started at 8 a.m. with initial conditions
as outlined in Chapter 6.
9.1.2 Undisturbed Soil

a. Metea sandy loam. The Michigan soil for which the
author was provided the most complete information and which
was modeled the most extensively was Metea sandy loam
(Figure 9.1)*. This soil exhibits the following

characteristics:

 

*Data compliments of Dr. Ray Kunze, Soil Science, Michigan

State University

117

118

 

Kl’=-‘)5C)OO cm.

 

4.0

 

3.0

‘l’ =-1I3 bar

“910‘?

20

 

KI’

Qt

 

 

90 .0'6 .1'2 .1'8 .2‘4 .36 .3'6 .4'2
Moisture Content, 9, m‘lm’

Figure 9.1. Moisture Characteristic Curve for
Metea Sandy Loam

119

8 at saturation = 37.5%

8 at 1/3 atm. = 16.0%
8 at 15 atm. = 8.5%
9H (immobile) = 2.0%

KS = 1.323-6 m/sec.

In addition, if the air-entry moisture content is taken as
90% of saturation, then

98 = 33.8%
and, from Figure 9.1, we = -.1405 m H20.

There was not enough data available to develop a mois-
ture adsorption curve. Hirschi and Moore (1980) stated
that, to obtain adsorption data, reasonable results are
derived by dividing the desorption data by 1.6. Using this
factor, simulations were based on we for adsorption a
-.1000 m H20.

The texture of Metea sandy loam was not available. The
author was advised to use the textural triangle, and the
following soil analysis was chosen:

55% sand
30% silt
15% clay

Further, a B horizon, starting at the 20-cm. depth, was

1

chosen as:
45% sand

30% silt
25% clay

1193

The bulk density for undisturbed, unstructured Metea
loam was set equal to 1650 Kg/mg, based on the saturation
moisture content of 37.5% provided in the data.

The surface roughness 20 was based on the values pre-
sented in Table 5.1, and was equated to 0.005 m for undis-
turbed soil.

The initial temperature at 8 a.m. was assumed a con-
stant value throughout the soil profile. An initial moisture
gradient was assumed for undisturbed soil. These initial
conditions were based on values in Table 9.1, Figure 9.1,
and several trial runs of the model.

b. Brookston clay loam. In order to study the '
effects of soil texture, a Michigan clay loam was also
modeled using the characteristic curve in Figure 9.2*. This

soil exhibits the following properties:

31.1% sand 8 at saturation a 33.6%
45.9% silt 8 at 1/3 bar = 22.1%
23.0% clay 8 at 15 bar = 13.0%
epwp = 11.3%

KS = 1.26E-5 m/sec

Bulk density at test = 1390 Kg/m3
At 90% of saturation, 8e = 30.2% and we = -.2690 m H20.
Dividing we by 1.6, wea = -.1680 m H20. The residual

 

*Data compliments of Dr. Ray Kunze, Soil Science, Michigan

State University

120

 

 

 

Q

In

C) i=~<£j,=-1EKDC“DICFU.

v:

Q

m i

'—-—-‘|’ =-1I3 bar
3.

 

I \
f3:

 

Q
,_

 

0 .06 £2 .18 .54 .50 .30
Moisture Content,e, rn’lm’

Figure 9.2. Moisture Characteristic Curve
for Brookston Clay Loam

121

(immobile) moisture content 8h was set equal to 11%. This
was based on a tested Opwp - 11.3% and values for clay loams
and silty clay loams recorded by Mualem (1978) ranging from
-9% to 28%.
The B1 horizon, starting at the 20-cm. depth, was
chosen as:
21% sand
44% silt
32% clay
Initial conditions and the surface roughness were chosen as
the same for the Metea loam.

Table 9.1. Soil conditions in northern locations
pertinent to model.

 

 

 

Soil Date T(°C) e(%) Depth(cm) Researchers
End
Nicollete of 18 S Cruse, et al.
loam May. (1980)
Spring 12.8 surface VanBavel and
to Hillel (1976)
24.0
Nicollete time 26 10
loam of 35 22 Allmaras, et al.
tillage 38 Below 22 (1977)
1965
time
of 22 20
tillage 26 30
1966
Sandy at 22 5 Ojeniyi and Dexter
loam F.C. (1979)
#1
Sandy at 12.6 S
loam F.C. 25.2 5
J;

 

 

 

 

 

 

 

122

9.1.3 Tilled Soil. Tillage operations during planting
substantially alter the structure of the soil. The model
already incorporates an expression for mean aggregate size
as a function of tillage passes (equation [65]) used in the
calculation of air permeability. But the greater effect on
the model is the reduction in the bulk density 08, which in
turn affects other key parameters. Bruce, et a1. (1977)
recorded an averageckss 1200 Kg/m3 in the top ten centi-
meters for a single pass of the tillage implement. Since
their Cecil sandy loam was similar to the Metea sandy loam,
this value was adopted for a single pass tilled zone in the
model. For two tillage passes through Metea, a DB = 1350
Kg/m3 was assumed.

Tillage also alters the moisture properties of the
soil. Hillel (1971) presented a qualitative discussion that
soil structure affects the w-8 relation, especially at low

suctions (see Figure 9.3).

    
 

compact

aggregated

 

Figure 9.3. Effect of Soil Structure on
Moisture Characteristics

123

For a given w near saturation, there is wide disparity
in 8. However, above we, and throughout most of the range,
the difference is generally diminished, and the two curves
converge at high suction. We can deduce from Ghosh (1980)
that the difference in we between the two curves in Figure
9.3 would be small. He argued that, if€u3= .90 , then eO-ee
is small, and we may be approximated as the intercept on the
line w=we (eleorB at 6/60 - 1. Therefore only one In; is
needed for a particular soil, independent of DB and.q3 (see

Figure 9.4).

€31€%,=1

---~————

    

veI—u—e-n-

 

Log (616.)

Figure 9.4. Ghosh Argument for One we

The we and wea for the undisturbed soil were assumed for all

tilled conditions as well.

124

The saturated hydraulic conductivity K5 is typically
given in the literature as a single value, even for studies
extending over variable bulk densities. Rawls, et al.
(1982) reported that KS, although dependent on ”B and other
factors, is much more dependent on texture, so investigators
typically ignore other factors, and record one KS for each
soil type.

However, to reflect a more accurate interdependence,
the 0B factor should be included.* Rawls, et al. cited an
expression of the form:

KS = f(sz)
Using this form, the following equation was used to calcu-

late “st of the tilled zone, knowing x5u of the undis-

turbed soil:

2650 - O
K = K Bt

- [86]
st su 2650 pBu

 

9.1.4 Meteorological Inputs. In order to simulate weather
conditions for central Lower Michigan, thirty-year statis-
tics were obtained from Bishop Airport in Flint, and are

shown in Tables 9.2 and 9.3.

 

*Dr. Ray Kunze of Michigan State University indicated till-
age could change £3 by a factor of ten, but very little

consistent data is presently available.

125

 

 

 

 

 

 

 

 

 

 

 

 

Table 9.2. Meteorological statistics from
Bishop Airport (1941-1970)
Temp. (°F) % Rel. Hum.
Month Daily Daily Monthly 0100 0700 1300 1900
Max. Min. Avg.
May 67.0 44.5 55.8 75 77 S4 56
June 77.0 54.6 65.8 79 80 57 58
Table 9.3. Diurnal meteorological statistics
for 1981 from Bishop Airport
Hour May Monthly Summary June Monthly Summary
Temp(°F) %Re1.Hum. Wind(mph) Temp(°F) %Re1.Hum. Wind(mph)
0100 49 73 6.7 61 79 6.4
0400 46 78 6.4 59 84 6.6
0700 48 76 7.0 62 80 8.2
1000 58 57 10.7 69 64 11.6
1300 62 50 11.8 74 57 12.5
1600 63 48 12.0 75 53 12.5
1900 60 52 9.2 72 56 9.8
2200 53 67 7.6 64 74 6.6

 

 

 

 

 

 

 

 

126

Global radiation data for the simulations were obtained
from Baker' and Haines (1969) and Baker and Klink (1975).
For simulating conditions in central Lower Michigan, statis-
. tics for East Lansing were gleaned from these reports and

are presented in Table 9.4.

Table 9.4. Mean weekly solar radiation over East Lansing,
langley/day.

 

 

 

Month Week 1953-1966______l2§2fi$212_522o
Avg. Clear sky 100% cloud Mean
May 1 490 725 265 487
2 469 740 150 465
3 504 680 200 501
4 544 700 260 539
June 1 550 735 160 544
2 532 760 250 529
3 565 750 240 555
4 580 760 235 552
5 548

 

 

 

 

 

 

 

 

Using the procedures outlined in Appendix D, meteoro-

logical inputs were calculated to simulate average May

weather conditions and hot, dry June conditions for Michi-

gan. These are shown in Table 9.5.

127

Table 9.5. Meteorological inputs used
in the simulations

 

 

 

 

 

 

Parameter Average Amplitude
Avg. May June Avg. May June
% Rel. 63 56 15 16
Humidity
Air temp 12.5 19.0 4.7 4.0
(°C)
Windspeed 4.1 4.1 1.3 1.3
(m/sec)
Global 50.0 50.0 685.0 774.0
Radiation
(W/mz)

 

 

 

 

The same diurnal cycle was imposed onto the model each day

for the duration of the simulation run.

2;; Simulation Results

9.2.1 Temperature and Moisture in Undisturbed Metea Sandy
Loam. Figure 9.5 shows the daily fluctuation in soil
temperature in undisturbed Metea loam. Successive noon
temperatures are cooling down, indicating sustained
meteorological inputs during simulation are insufficient to
maintain the temperature level dictated by the 8 a.m. start
condition. Temperatures seem to be approaching equilibrium
by the third day. Figure 9.6 exhibits the corresponding
moisture conditions. The entire profile is reduced in
moisture content, producing a steeper gradient by the third

day. Although moisture depletion in the sublayers is

Depth , cm.

128

 

 

 

 

 

 

 

 

 

Temperature ,°C
06 a 10 1_2 1.4 Is 18
F‘“\ 1
N 021;Ҥ( Tl /
C2». \\
‘9’5z\ I
I0 i
l
9.-..-- l
I
I
l
t)
19
ml
0 21
N (I I
00 l' ‘
I; \
ll ‘\
I
I}! .
:I \
o J! \
m

 

Figure 9.5.

Bulk density=1650 kg I m3

Temperatures in Undisturbed Metea Sandy
Loam for Average Michigan May Climate

129

Moisture Content , e , m’I m’
00 A .04 . .08 . .12 J .10

 

. -2.o . 2.4 .
\\ \
Ni \\ \
\
I \
ID \ I

 

 

 

 

‘ \
S2 i
\ \ Days
. ‘ 1
2
g ‘I3 \
f.‘ ‘i \
8- I
o I |
O i , m
N I ‘
i
‘I |
‘. |
'. I
o - ‘. l
m

 

Bulk density =1650 kglm3

Figure 9.6.

Noon Moistures in Undisturbed Metea Sandy
Loam for Average Michigan May Climate

130

slowing down, drying of the top layers is accelerating. The
drying t0players have less thermal inertia, thereby cooling

down further at night.

9.2.2 Temperature and Moisture in Undisturbed Brookston
Clay Loam. Figures 9.7 and 9.8, showing temperature and
moisture profiles in Brookston clay loam, can be compared to
Figures 9.5 and 9.6, the counterparts for Metea sandy loam.
The sandy loam has a higher thermal conductivity than the
clay loam, but its higher bulk density creates a higher
thermal inertia with fewer air pockets. Thus, the two
appear to cancel each other in the temperature regime.

The clay soil retains more moisture in the top layers
during the first day, but at the expense of deeper layers.
By the third day, although the surface in the clay soil has
dried out further, moistures below the surface are
approximately the same for both soils. This surface
condition probably can be attributed to a larger vapor flow
component present near the surface for the lower bulk
density.

The B1 horizon of the clay soil restricted moisture
movement, as evidenced by the moisture plot distortions in
Figure 9.8. The evaporative demand thus drew more moisture
from the 10- to 20-cm depth in the clay soil than from the
sandy soil since moisture at deeper layers was not as

available here as it was in the sandy soil.

131

Temperature ,°C

6 8 10 12 1,4 16 Is 20

 

 

 

 

 

 

 

 

 

0" L} |
N
N 05 i
\\ I
In 9; . t ,1
0. / 0
<9) / 0°
<"\ , e
l /
$2 --___ “__._ _ .
. Days
6 | ’ 1
., ‘gI x/’ 2———
5 17.. I/ 3"""
8 2|
0 <1 / \
o 00" i
CU l : ‘
I: l
I! I
' i
II
o .
00

 

Bulk density =1390 kg/m’

Figure 9.7. Temperatures in Undisturbed Brookston Clay
Loam for May

Depth , cm.

132

Moisture Content , e, m’l m3

 

 

 

 

 

 

 

 

.04 .O§ .12 16 . 29 . 2f) . 25
o ._._\ / -
N \‘ \
\ \ \
\
\ \
ID : ~—
\\ \ \
I‘ \
$2 ! i
I Days
I ‘ 1
l l2
I
I3 I
I
i l
8 L l
\\\ \
\
\\\\
O
F)

 

 

Bulk density =1390 kg/m3

Figure 9.8. Noon Moistures in Undisturbed

Brookston Clay Loam for May

133

9.2.3 Comparisons Between Undisturbed and Tilled Metea
Loam. Figure 9.9 compares moistures at several depths for
undisturbed Metea loam and after a single tillage pass
creating a lO-cm. surface layer at 1200 Kg/m?. The diurnal
moisture cycle is obvious with undisturbed soil, showing
fluctuations highest at the surface and decreasing with
depth. Peak-to-peak fluctuations these first three days
ranged from 2% to 6% at the surface, and from 1% to 2% at
five centimeters. Bruce, et al. (1977) reported identical
fluctuations at these depths for their sandy loam. Hillel
(1977), in his computer model of a fine sandy loam,
exhibited a 10% fluctuation at the five millimeter depth.

Tilling the soil to a bulk density of 1200 Kgim3'
appears to significantly alter these diurnal cycles. The
surface moisture dries quickly to virtually zero (artifi-
cially held to 1% by the model), and increases only 2 to 3%
during the early morning hours. Meanwhile, moisture at the
5-cm depth steadily increases over the three-day period.
This phenomenon is explained in Figure 9.10.

Figure 9.10 shows the three-day history of the noon
moisture profile. The moisture discontinuity created at the
interface at 8 a.m. of the first day readjusts itself to a
continuous gradient by the third day. The t0p layers of
soil dry out, but the rate of drying slows by the third day.
Moistures at depths below 3 cm actually increase over the 3-
day period. Moisture collects below the interface as verti-

cal movement from below proceeds but movement across the

134

wocmm mmumz How magnum mQOme> um mmuoumfloz

pwaawa pom omnusumflpco mcflummeoo Smog

 

  

 

.m.m mesons
9.301
8.6.1. n.8,», . fat... no
\I, :11), .1---)

x xi: \ I t .

uumtsm. II /I\ I \ .
I.\ '.0
m xmo N came :30 .7

159.0 me 1

00923
.E0 0

.Eu 9

.cCuvmu\

 

    

MCE? o

l
‘
“'I||
\‘ll‘

 

:80 99 00.5.

9 68336:: ill

r"'

 

 

.
_
_
_
_
_
.1
_
—
_
.
_
_

80'

,w [,w ' e 'iuaiuog eJnisgow

er T éI‘ '

I

03'

 

135

Moisture Content ,6 , m’/ m’
0 .04 .08 , .1g 1 .10 _20 . 2.4 gala

 

 

 

 

 

 

O _________ \
N _____—_-_-_.::::.2‘ 333
tilled - 1200 kglm'
:- \\\</// undisturbed -1650 kglma
6
U
1:.
4.:
CL
0)
o
o
o
00

 

Noon Moisture Profiles in Metea Loam with

Figure 9.10.
a Tilled Bulk Density of 1200 Kg/m

136

interface is impeded. Moisture content has reached a uni-

form level throughout the profile below ten centimeters.

Figure 9.11 compares temperatures between the undis-
turbed and tilled condition. Tilling warms up the entire
profile by the third day. The lO-cm. tilled layer warms up
faster, creating a slight discontinuity in the temperature

gradient at the lO-cm depth.

9.2.4 Tilled Density and Tillage Passes Comparisons in Metea
Loam.

Figure 9.12 compares the effects of bulk density. A
bulk density of 1200 Kg/m3 is representative of a single
pass of a tillage implement, whereas 1350 Kg/m3 is more
representative of two passes. Here again moisture trends
reverse at a depth of about 3 centimeters. Below this point
moisture contents by the third day are higher for the lower
bulk density, possibly due to lower liquid movement. Above
this point moisture is lower for the lower bulk density,
probably due to higher vapor flow, the predominant mechanism

in dried surface layers.

9.2.5 Minimum Tillage Furrow Shape Comparisons. Figure
9.13 shows several furrow dimensions and shapes which were
modeled using the grid of Figure 6.3. The oval at the base
of each furrow is the planted seed. Different combinations
of the small elements in the upper right hand corner of the

grid in Figure 6.3 were used to simulate the furrows.

Depth , cm.

137

Temperature °C

 

 

 

 

 

 

 

 

6 10 22 26 39
O k ‘ { ,1
00‘ ,I/
\(9 \ /
,0 In \ R /
C) ,{L \‘/,r
tilled
Q ///)\\\\
“ undisturbed
I I
I I
I I
I I
’ I
i I
O . e
N 'T 1'
g ' ---tilled
, ' 1200kg/m’
' : —-undisturbed
: | 16:30kg/m3
I
‘ I
O L L
m

 

Figure 9.11.

Third Day Temperatures in Metea Loam
Comparing Undisturbed and Tilled
Condition

138

Noon Moisture Content, 9, m’lm’

 

 

 

 

 

 

0 .04 . .08 . .12 1 .19 .29 _25 .218
01"?" ———————— ‘ I
\—-—_ \
N- ““““““ 1
\\
G) \\
p
‘n z
(D
E?
tilled 3,1
9 ///A\V
E- undisturbed
k)
c.
4.:
o.
0
O
O
N Bulk density kg/m’
1200—fl
1350 , I
I I
I I
l I
o I I
m

 

 

Figure 9.12.

Comparing Tilled Bulk Densities

Third Day Moisture Profiles in Metea Loam

139

The results of the simulations of these furrow shapes
are exhibited in Figure 9.14. Moisture content at the seed
point decreases at the same rate for all furrow shapes for
the undisturbed condition. However, the seed point stays
wetter by 1 to 2% when at the base of a tilled furrow as
compared to when in the undisturbed profile. This, again,
is due to moisture movement from lower layers up to the
seed, but impedance to moisture movement through the furrow.
Although there is a difference of up to 1% moisture content
between the wide rectangular and narrow furrows, there is
essentially no difference between the two narrow) profiles.
This might suggest that moisture at the seed is more

dependent on furrow width than furrow shape.

 

I;

1200 1200 1200 §1200 i

//K,\\

 

 

/VK3x.
Aocwx

 

I<——UI——)

 

 

 

 

 

 

 

50
l-e—B {Ala-I |<—4 —>I I<-4-a-I I64 ->I
- a. b. c. d.
wide narrow narrow narrow V
rectangle rectangle V modified

all units = cm.

Figure 9.13. Minimum tillage furrows modeled.

Figure 9.l3(d) shows a narrow V furrow simulating a

condition where some loose soil fell into the furrow in

140

0umewau >0:
0m0u0>< new 8004 hocmm 00002 suds m0aamoum souusm
00HHNB H0u0>0m How A80 my unwom 000m um 0usumwoz .va.m 0usmfim

 

 

 

 

 

 

9501
VN we m.VN.0%.m.VN.©.F.mO
m >00 N >00 :30 .

. w

.0
v m.
1'
s n.
..O a
8 O
o
u “u
030.6 r. «e.
p0ne3mfica/ .8 W

>>>otmc.u .

03:309. .. .9

>>otmcfi .L
03:300.. \ .9 m...
023.0 . w...

6 m6

. . .O

M r

141

front of the seed firming wheel and was firmed immediately
over the seed to 1500 Kg/m3, followed by loose soil back-
fill. This furrow condition was simulated on the computer
and compared with the narrow V shape of Figure 9.13(c).
Three nodes of the grid, all on the furrow centerline, were
compared. The node at the surface for both conditions dried
to 1% and remained there for the 3-day period. The node
representing the seed did not vary more than 0.2% between
conditions for the simulation period.

The largest difference was experienced with the node at
the interface of the two densities. In the modified furrow,
this node, in the first day, was 4 to 5% higher in moisture
content than the narrow V furrow. However, by the end of
the third day, this had reversed where the modified furrow
was 3 to 4% lower than the narrow V.

Apparently the interface created a moisture barrier
early in the period, and water accumulated at the interface.
Later, after continuity had been restored, the 1500 Kg/m3
section of soil dried out further than had it been at 1200
Kg/mi3. The narrow V of Figure 9.13(c) had built up a
barrier to liquid flow earlier in the period, where as the
modified. condition allowed additional moisture to be
”trapped“ in the triangular section to be lost later. The
two conditions simulated had no effect on moisture at the

seed, suggesting it had received adequate moisture from

lower depths throughout the 3-day period.

142

9.2.6 Varying Meteorological Conditions. Figures 9.15 and
9.16 show the effect different meteorological conditions
have on the temperature and moisture profiles. Average May
and hot June conditions for Michigan, as outlined in Table
9.5, were used in this simulation. With average air temper-
atures increased by 6.5°C, the soil temperature by the third
day has increased by 6°C throughout most of the profile.
However, with a corresponding reduction of average air rela-
tive humidity by 7%, most of the moisture profile remains

unchanged. Only the top one centimeter has dried further.

9.3 Summary.
This chapter exhibited the results of simulations of

soil conditions comparing soil types, structure, tilled
furrow shapes and dimensions, and meteorological conditions.
Chapter 10 presents concluding remarks and recommendations

for future work.

143

Temperature , °C
5 10 1g 29 2;) 30 3;:

 

 

 

 

 

 

 

 

 

0
(Vs
In
0 ,/ tilled
v- \/ ///>\\\\
. ’l undisturbed
E / I
U. / I
.c: I I
E
0 / l
0 I I
o J I _
N I I I
average / I hot June
May / l \I\ days
days I I
I I
I I
I
O I
m

Figure 9.15. Third Day Temperature Prgfiles in Metea
Loam Tilled to 1200 Kg/m Exposed to
Varying Meteorological Conditions

144

Noon Moisture Content ,9 , m’/ m3
0 .04 .08 .12 .16 .20 .24 .28

 

 

 

 

 

 

 

 

O I a l 1‘ A 4 j l 1 1 L 1 4 n
LLé—z—L- ::::: - --/L’ hOt June days
Nq average May days/ “r
In
(n
>
, Z
O tilled ~1200 kglm’
'- f///>\\\\
. undisturbed -1650 kglm’
E
(J
1:
.p
o.
o
o
O
N
I
I
I
l
I
I
O l
m

Figure 9.16. Third Day Moisture Profiles in Metea Loam
Tilled to 1200 Kg/m3 Exposed to Varying
Meteorological Conditions

10. CONCLUSIONS AND FUTURE WORK
lggl Summary

This study produced the following, in reference to the
seed-soil system model of Figure 1.1:

a. the two-dimensional soil physical model

b. the meteorological model

c. the soil surface interface tying the boundary con-
ditions of the soil model to the meteorological inputs

The model was based on the Galerkin approximation of
the finite element method of numerical solutions to simul-
taneous differential equations relating moisture and temper-
ature interactions. The computer program was developed in
the FORTRAN language and run on the IBM 360 computer at
Delta College. It was validated with actual field data
extracted from the literature, and with few modifications,
exhibited good agreement with the data.

Although simultaneous interactions of two phenomena and
time-dependent problems have previously been modeled with
the finite element method, the model of this study is unique
due to the combination outlined below:

1. Uses Galerkin approximation to finite element numerical

solution
. Is two-dimensional and time—dependent
. Includes simultaneous moisture and temperature interactions
Includes heterogeneous soil grid

. Interfaces with real meteorological model

0‘ 0" 9 N N
e

Applies to seed furrow and tillage practice studies

145

146

Simulation results shown in Chapter 9 demonstrate that

the model can be used to study the following parameters and

their effects on soil moisture and temperature:

01 uh I» N 5..-
0

Soil texture of AP and B1 horizons

Soil structure and tillage effects on structure
Furrow geometry

Meteorological conditions

Albedo

10.2 Conclusions

The objective of this study was to develop a realistic

model of soil moisture and temperature, not to develop new

theory.

Nevertheless, several conclusions can be drawn from

the experience gained in developing the model.

Moisture diffusivity due to vapor flow induced by
moisture gradients cannot be ignored, and in fact
is a significant factor at low moisture content in
the surface layers, especially during mid
afternoon.

The "mass enhancement factor" added to the soil
thermal conductivity term in equation [57] fully
accounts for heat convection through the soil due
to vertical air movement. A separate air convec-
tion term using the vertical air velocity gradient
3V/82 in addition to the mass enhancement factor
produces excessive heat flow conditions.

The soil air mass flow term due to vertical air

density gradients, obtained froml equation [68],

147

produces excessive moisture movement. Whereas the
theory underpredicts this phenomenon, the model
developed in this study overpredicts.

d. The soil moisture diffusivity due to temperature
gradients is several orders of magnitude smaller
than that due to moisture gradients. Only at very
dry conditions below the wilting point does De sink

to a level where DT becomes significant.

lg;; Future £955

The models developed from this work could be considered
a good start towards the attainment of the long-range goal
identified in Section 1.3. Yet a great deal of work needs
to be done before the seed-soil system model of Figure 1.1
can be realized. This work can be categorized as follows:

1. Second validation of the present models

2. Refinement of parameters of present models

3. Expansion of utility of present models

4. Development of three-dimensional seed model

5. Development of three-dimensional soil model with

moisture, temperature and oxygen parameters.

6. Linking of all models into seed-soil system model.

10.3.1 Second Validation

Although the soil model has shown good agreement with
field data from one study, a second validation with a dif-
ferent soil under different weather conditions would add

credibility to the model. Field data for such a validation

 

148

would need to supply the inputs listed in Sections 4.7 and
5.3. At the time of this study, sufficient published data
were not available for a second validation. The field study
by Jackson, et al. (1973, 1975) came the closest to supply-

ing all inputs, but several key parameters were missing.

10.3.2 Refinement of Parameters

Chapter 8 dealt with model refinements conducted during
this study. Although refining a model is virtually a con-
tinuous process as experience with it accumulates, several
specific suggestions are listed here which could improve the
accuracy of the model.

a. Surface conditions. Section 8.2.1 explained the
physical constraint imposed on the surface temperature in
order to suppress oscillations in the finite element solu-
tion. The surface condition requires further development
and may justify a separate model to be interfaced with the
soil model.

Closely coupled with this are the evaporation and sen-
sible heat flux terms in the surface energy balance equation
explained in Section 8.2.3. Although these were taken
directly from the literature, closer study of their inter-
actions would possibly benefit the overall seed-soil system
model.

b. Soil moisture diffusivity De. The moisture diffu-
sivity De was used with relationships based on the most

advanced theory to date. However, these relationships,

149

embodied in equations [44] and [47], are still only approxi-
mations, with limitations to their applicability to all soil
classes. The expression for saturated hydraulic conducti—
vity as a function of density, equation [86], is also an
approximation. The soil model will have to incorporate
improvements in these relationships as theory is refined in
the future.

c. Hysteresis. Although the model incorporates hys—
teretic effects in the soil, virtually no field data are
available to verify this parameter. The development of
adsorption data in the field for modeled soils is necessary

to quantify this part of the model.

10.3.3 Model Expansion

a. The use of the soil model to simulate conservation
tillage or minimum tillage practices will be enhanced with
the incorporation of a surface model which represents corn
stalks and other surface residue. It will affect such para-
meters as albedo and surface roughness, and would create a
separate layer of finite thickness with a unique moisture
diffusivity, thermal conductivity and heat capacitance.
Literature on this topic is still sketchy, and considerable
research is implied.

b. The model in its present form accepts a diurnal
Cycle of weather conditions which is the same every day.
More realistic conditions can be simulated if the average

and amplitude of each cycle can be changed from day to day.

150
c. Presently the model does not accept a rainfall.
The incorporation of this weather condition is important
when expanding the model's use to the study of seed response

to oxygen deficiency.

10.3.4 The Three-Dimensional Seed Model.

Referring back to Figure 1.1, a three-dimensional seed
biophysical model must be developed to determine the hourly
viability of the seed in response to its environment. This
model must represent the biophysical processes within a
given seed which would include temperature, moisture, and
oxygen cause-effect relationships. Throughout the past few
decades, a number of investigators have studied the biophy-
sical processes of germinating seeds. The author has liter-
ature on seed studies dating to 1916, and recent partial
models advanced, exemplified by that by Waggoner and
Parlange (1976), would indicate at least a crude model could
be developed, drawing from submodels found in the litera-

ture.

10.3.5 The Three-Dimensional Soil Model.

A spherical soil shell wrapped around the spherical
seed model needs to be developed to interact directly with
the seed. According to Figure 1.1, at each time step, this
soil shell receives its moisture and temperature distribu-
tion from the two-dimensional soil model. The shell's

radius extends out from the seed center only to the extent

151

of the seed's influence. This model would include a spheri-
cal seed-soil interface providing impedance to heat, mois-
ture and oxygen exchange between the seed and soil.

The author is aware of several studies of seed imbibi-
tion as affected by seed-soil contact area. Earlier inves-
tigations are represented by Collis-George and Hector (1966)
and Hadas (1969). More recent work is represented by

Bruckler (1983) and Boiffin, et a1. (1983).

10.3.6 The Seed-Soil System Model.
Once the various component models have been verified
and refined, a final linking together is necessary to obtain

a harmonious flow as outlined in Figure 1.1.

LIST OF REFERENCES

LIST OF REFERENCES

1. Allmaras, R.R. et al.
1977. Surface energy balance and soil thermal property
modifications by tillage-induced soil structure.
Minnesota Agric. Exp. Sta. Technical Bulletin

306.

2. Arkin, G.F. et al.

1980. Forecasting grain sorghum yields using simulated
weather data and updating techniques.
Transactions of American Society of Agricultural
Engineers, 23(3):676-680.

3. Baker, D.G. and D.A. Haines.

1969. Solar radiation and sunshine duration
relationships in the North-Central Region and
Alaska. North Central Regional Research
Publication 195, Technical Bulletin 262. Agric.
Exp. Sta. Univ. of Minnesota.

4. Baker, D.G. and J.C. Klink.
1975. Solar radiation reception, probabilities, and
areal distribution in the North-Central Region.
North Central Regional Research Publication 225,
Technical Bulletin 300. Agric. Exp. Sta. Univ.
of Minnesota.

5. Boiffin, J. et al.

1983. Role des proprietes physiques du lit de semences
sur l'imbibition et al germination. III. Valeur
previsionnelle d'un modele d'imbibition au champ
et caracterisation des lits de semences.
Agronomie 3(4):291-302.

6. Bresler, E. et al.

1969. Infiltration, redistribution and subsequent
evaporation of water from soil as affected by
wetting rate and hysteresis. Soil Science
Society of America Proceedings, 33:832-839.

7. Bruce, R.R.
1972. Hydraulic conductivity evaluation of the soil
profile from soil water retention relations.
Soil Science Society of America Proceedings,

36:555-561.

8. Bruce, R.R. et al.
1977. Diurnal soil water regime in the tilled plow
layer of a warm humid climate. Soil Science
Society of America Journal, 41:455-460.

152

153

9. Bruckler, L.

1983. Role des proprietes physiques du lit de semences
sur l'imbibition et la germination. I.
Elaboration d'un modele du systeme terre-graine.
II. Controle experimental d'un modele
d'imbibition des semences et possibilites
d'applications. Agronomie 3(3):213-222 and 223-

232.

10. Cary, J.W. and S.A. Taylor.
1962a. The interaction of the simultaneous diffusions of
heat and water vapor. Soil Science Society of
America Proceedings, 26:413—416.

11. Cary, J.W. and S.A. Taylor.
1962b. Thermally driven liquid and vapor phase transfer
of water and energy in soil. Soil Science
Society of America Proceedings, 26:417-420.

12. Cary, J.H.
1965. Water flux in moist soil: thermal versus suction

gradients. Soil Science, 100 (3):l68-175.

13. Collis-George, N. and J.B. Hector.

1966. Germination of seeds as influenced by matrix
potential and by area of contact between seed and
soil water. Australian Journal of Soil Research,

4:145-164.

14. Cruse, R.M. et al.
1980. A model to predict tillage effects on soil
temperature. Soil Science Society of America

Journal, 44:378-383.

15. Cushman, J.H. and Kirkham, D.
1978. A two-dimensional linearized view of one-
dimensional unsaturated-saturated flow. Water
Resources Research, 14(2):319-323.

16. Cushman, J.H. et a1. .

1979. A Galerkin in time, linearized finite element
model of two-dimensional unsaturated porous media
drainage. Soil Science Society of America
Journal, 43:638-641.

17. DeVries, D.A.
1958. Simultaneous transfer of heat and moisture in
porous media. Transactions of American
Geophysical Union, 39(5):909-915.

154

18. DeVries, D.A.

1963. Thermal properties of soils. In: VanWijk,
Physics of Plant Environment. North-Holland
Publishing Co. Amsterdam, Div. of John Wiley &
Son, NY:210-235.

19. Farrell, D.A. et a1.
1966. Vapor transfer in soil due to air turbulence.

Soil Science, 102:305-313.

20. Fukuda, H.
1955. Air and vapor movement in soil due to wind
gustiness. Soil Science, 79:249-258.

21. Gebhart, B.
1961. Heat Transfer. McGraw-Hill Book Co., Inc. New

York.

22. Ghosh, R.R.
1980. Estimation of soil-moisture characteristics from
mechanical properties of soils. Soil Science,

130(2):60-63.

23. Gray, W.G. and G.F. Pinder.
1974. Galerkin approximation of the time derivative in
the finite element analysis of groundwater flow.
Water Resources Research, 10(4):821-828.

24. Gustafson, R.J. et a1.
1979. Temperature and stress analysis of corn kernel-
finite element analysis. Transactions of
American Society of Agricultural Engineers,
22(4):955-960.

25. Hadas, A.
1968. Simultaneous flow of water and heat under
periodic heat fluctuations. Soil Science Society
of America Proceedings, 32(3):297-301.

26. Hadas, A.

19693. A comparison between predicted and measured
values of the thermal conductivity of a moist
soil under steady and fluctuating thermal
regimes. Israel Journal of Agricultural
Research, 19(4):151-159.

27. Hadas, A.
1969b. Effects of soil moisture stress on seed
germination. Agronomy Journal, 61(2):325-327.

155

28. Hadas, A.
1977a. Evaluation of theoretically predicted thermal
conductivities of soils under field and
laboratory conditions. Soil Science Society of
America Journal, 41:460-466.

29. Hadas, A.
1977b. Heat transfer in dry aggregated soil: I. Heat
conduction. Soil Science Society of America
Journal, 41:1055-1059.

30. Haghighi, K. and L.J. Segerlind.
1978. Computer simulation of the stress cracking of
soybeans, American Society of Agricultural
Engineers paper no. 78-3560. St. Joseph, MI.

31. Haghighi, K.
1979. Finite element formulation of the thermo-hydro
stress problem in soybeans, Unpublished Ph.D.
Thesis, Michigan State University, East Lansing,
MI.

32. Hanks, R.J. and N.P. Woodruff.
1958. Influence of wind on water vapor transfer through
.soil, gravel, and straw mulches. Soil Science,
86:160-165.

33. Harper, L.A. et al.
1976. Soil and microclimate effects on Trifluralin
volatilization. Journal of Environmental
Quality, 5(3), July-Sept.

34. Hillel, D.
1971. Soil and Water: Physical Principles and
Processes. Academic Press, New York.

35. Hillel, D.
1977. Computer Simulation of Soil-Water Dynamics.
International Development Research Center,
Ottawa, Canada. -082e:35-78.

36. Hirschi, M.C. and I.D. Moore.
1980. Estimating soil hydraulic properties from soil
texture. American Society of Agricultural
Engineers paper no. 80-2523. St. Joseph, MI.

37. Hollman, J.P.
1976. Heat Transfer. McGraw-Hill Publishing, New York.

38. Jackson, R.D. et al.
1973. Diurnal soil-water evaporation:time-depth-flux
patterns. Soil Science Society of America
Proceedings, 37:505-509.

156

39. Jackson, R.D. et a1.
1975. Heat and water transfer in a natural soil
environment. In: Heat and Mass Transfer in the
Biosphere. John Wiley & Sons, New York:67-76.

40. Judah, O.M. et al.
1975. Finite element simulation of flood hydrographs.
Transactions of American Society of Agricultural
Engineers, 18(3):518-522.

41. Jury, W.A. and J. Letey.
1979. Water vapor movement in soil: Reconciliation of
theory and experiment. Soil Science Society of
America Journal, 43:823-827.

42. Kimball, B.A. et al.
1976. Comparison of field-measured and calculated soil-
heat fluxes. Soil Science Society of America
Proceedings, 40:18-24.

43. Laroussi, C. et al.
1975. Experimental investigation of the diffusivity
coefficient. Soil Science, 120 (4):249-255.

44. Luikov, A.V.
1966. Heat and Mass Transfer in Capillary-Porous
Bodies. Pergamon Press Ltd. London.

45. Malik, R.S. et al.
1979. Physical components of the diffusivity
coefficient. Soil Science Society of American
Journal, 43:633-637.

46. Mark's Mechanical Engineers' Handbook. McGraw-Hill Book
Co., New York.

47. Merva, G.
1975. Physioengineering Principles. AVI Publishing
Co., Westport, Conn.

48. Misra, R.N. and J.H. Young.

1979. The finite element approach for solution of
transient heat transfer in a sphere.
Transactions of American Society of Agricultural
Engineers, 22(4):944-949.

49. Misra, R.N. and J.H. Young.

1980. Numerical solution of simultaneous moisture
diffusion and shrinkage during soybean drying.
Transactions of American Society of Agricultural
Engineers, 23(5):1277-1282.

157

50. Misra, R.N. et a1.

1981. Finite element procedures for estimating
shrinkage stresses during soybean drying.
Transactions of American Society of Agricultural
Engineers, 24(3):751-755.

51. Mualem, Y.
1978. Hydraulic conductivity of unsaturated porous
media: Generalized macrosc0pic approach. Water
Resources Research, 14(2):325-334.

52. Ojeniyi, 8.0. and A.R. Dexter.
1979. Soil structure changes during multiple pass
tillage. Transactions of American Society of
Agricultural Engineers, 22(5):1068-1072.

53. Pall, R. and Mohsenin, N.N.

1980. A soil air pycnometer for determination of
porosity and particle density. Transactions of
American Society of Agricultural Engineers,
23(3):735-741, 745.

54. Pall, R. et a1.

1980. Comparison of one-dimensional flow simulation by
two numerical methods. American Society of
Agricultural Engineers paper no. 80-2526, St.
Joseph, MI.

. 55. Philip, J.R. and D.A. DeVries.

1957. Moisture movement in porous materials under
temperature gradients. Transactions of American
Geophysical Union, 38(2):222-231.

56. Potter, M.C. and J.F. Foss.
1975. Fluid Mechanics. John Wiley & Sons, New York.

57. Rawls, W.J. et al.
1982. Estimation of soil water properties.
Transactions of American Society of Agricultural
Engineers, 25(5):1316-1320, 1328.

58. Rose, C.W.
1968a. Water transport in soil with a daily temperature
wave: I: Theory and experiment. Australian
Journal of Soil Research, 6:31-44.

59. Rose, C.W.
1986b. Water transport in soil with a daily temperature
wave: 11. Analysis. Australian Journal of Soil

Research, 6:45-57.

60. Segerlind, L.J.
1976. Applied Finite Element Analysis. John Wiley &
Sons, New York.

158

61. Sepaskhah, A.R. and L. Boersma.
1979. Thermal conductivity of soils as a function of
temperature and water content. Soil Science
Society of America Journal, 43:439-444.

62. Townsend, J.S. and P. Migchels.
1981. Albedo of wheatfields under five tillage systems
in Manitoba. American Society of Agricultural
Engineers paper no. 81-1014. St. Joseph, MI.

63. Tscheschke, P.D. and J.R. Gilley.
1979. Status and verification of Nebraska's corn growth
model - CORNGRO. Transactions of American
Society of Agricultural Engineers, 22(6):1329-
1337.

64. Van Bavel, C.H.M. and D.I. Hillel.
1976. Calculating potential and actual evaporation from
a bare soil surface by simulation of concurrent
flow of water and heat. Agricultural
Meteorology, 17:453-476.

65. VanWijk, W.R.
1963. Physics of Plant Environment. North-Holland
Publishing Co., Amsterdam, Division of John Wiley
& Sons, New York.

66. Waggoner, P.E. and J. Parlange.
1976. Water uptake and water diffusivity of seeds.
Plant Physiology, 57:153-156.

67. Woodside, W. and J.M. Kuzmak.
1958. Effect of temperature distribution on moisture
flow in porous materials. Transactions of
American Geophysical Union, 39(4):676-680.

68. zienkiewicz, O.C.
1971. The Finite Element Method in Engineering Science.
McGraw-Hill, London.

APPENDICES

APPEND I X A
COMPUTER FLOW D I AGRAM
The flow of information, processing and key decision
points of the computer program are outlined in diagram form.

This serves as a ”bird's eye" view of the total program.

 

v

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

READ DO

Soil and Air on each level ‘1
Properties, ‘
Functions \L

i CALCULATE

air velocity

CALCULATE gradient BV/az
parameters, -
matrices

¢ END DO LOOP

 

 

 

 

I W < __IL

on each node

* DO ,
\L on each element

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mu V
Node locations, CALCULATE
depths, properties,3V/az,
initial conditions evaporation and heat
J! at surface
(END DO Loop]L \L
’ FORM

 

element conductance
matrix [ke] and

 

 

 

 

 

 

    
 

   

 

DO A A heat vector {f9}
on each element for equation [42]
I
READ boundary Yes

   

Element grid, -1eme
boundary elements,

boundary conditions

  

 

 

 

ADD
\II heat f low‘ term
tolife} I

 

 

 

 

 

 

 

 

 

 

( END no LOOP}

DO

on each time step (B)
INITIALIZE

Global Matrix [K] and

Global Vector {F}
[

 

 

 

 

 

 

 

 

Figure A.l Computer Flow Diagram

__E’? °

 

INSERT
[ke] into Global
Matrix [K] and
{re} into Global
Vector {F

W’

END DO L00P“‘
C r

 

 

 

 

 

 

 

SOLVE
Decompose [K] and
{F} and solve for T

I

INITiALIzE
[K] and {F}

I

DO 12v,

 

 

 

 

 

 

 

 

 

 

on each element

I

CALCULATE
Properties

I

 

 

 

 

 

FORM
element diffusivity
matrix [ke] and
eva oration vector
{f3 for equation [43]

 

 

 

 

 

 

 

 

ADD
evaporation
term to {re}

 

 

 

 

INSERT

 

{fe}into {F}

[ke] into [K] and

 

 

 

 

END DO LOOP#/i

 

I

 

 

SOLVE

Decompose [K] and

F} and solve
for e

 

 

 

I

DO
on each node

I

%

 

 

PRINT
all nodal values

 

 

 

 

 

C END Do LOOP}———

 

 

 

l\
(TEND DO LOOP#/

[END

 

APPENDIX B
COMPUTER PROGRAM LISTING
The main program and its eight subroutines are listed,
complete with sufficient documentation on variables and data

entry for the reader to engage in future work.

MAIN PROGRAM DOCUMENTATION

Computer Variable

Explanation

 

B and CI vectors

BC

ELK
BULK

DIEM

DIET

vector used to arrange global field
variables in column form conducive to

solutions in subroutines DCMPBD and

SLVBD

2 1 1
1 2 1
1 1 2

amplitude of horizontal air pressure
waves at surface, sec"l

area of element considered

elements in the BC matrix

common matrix used to build element con-
ductance and diffusion matrices

bulk density of tilled soil, kg/m3

soil bulk density, kg/m3, generally of
undisturbed soil

heat capacity of total soil volume and
air-vapor mixture, product of subroutine
CAP

element moisture diffusivity 09' product
of subroutine DIFFM

element moisture diffusivity DT, product

of subroutine DIFFT

DVDDZ

DVDZ

DVS

ECM

EDM

EEV
EF

ELEV

GCM, GHV

GGF (no. nodes),
GGSM (no. nodes ,
bandwidth)

vector of BVD/az values for all depths
in grid

aV/az value for element considered,
product of subroutine VEL

volume surface mean diameter of soil
aggregates of element considered
surface evaporation vector

element heat conductance matrix [8].
3x3

element moisture diffusion matrix [R].
3x3

surface element evaporation force vector
intermediate force vector for surface
elements used to build final force

vectors EHV and EEV
build final force vectors EHV and EEV

surface element heat flux force vector
elevation, meters, where wind, vapor
pressure, etc. were measured

frequency of horizontal air pressure
waves at surface, sec-4'

global conductance matrix and global
force vector used to store variables in
preparation for A vector

global stiffness matrix and global force
vector used to solve equations in sub-

routines DCMPBD and SLVBD

GPHI

- HCON

HCONS

HONS

IS?

IT

ITIMET
JEND

JG?

JGSM

field variables (temperature or moisture
content) produced from solutions‘in
subroutines DCMPBD and SLVBD

element hydraulic conductivity, product
of subroutine DIFFM

saturated hydraulic conductivity for
element considered, m/sec

saturated hydraulic conductivity of
tilled soil, m/sec

index to indicate status of soil
properties:0 - all elements same, 1 -
some elements are different

time increment between successive steps,
minutes

total time of simulation run, minutes
a pointer indicating the last storage
location for [K] in {A}

a pointer indicating the last storage
location for {¢} in [A]

a pointer indicating the last storage
location for {F} and [A]

total no. of horizontal levels in grid
(surface included)

bandwidth ((R+1)NDOF) used in
subroutines DCMPBD and SLVBD (see
Segerlind (1976))

total no. of elements

NI, NJ, NR
NIP

NLEVl
NLEVZ

NS vector

PSI

PSIE, PSIEA

QAV. 0AM

no. of elements which have boundary

conditions

.no. elements having soil structure

different from the majority. Program
allows two different sets of properties
(i.e. tilled vs. undisturbed)

no. elements with texture different from
majority (usually Bl horizon)

three nodes of element considered

no. implement passes

first node at specific depth considered
second node at specific depth considered
total no. of nodes

no. nodes having initial 8 or T
different from majority

no. global degrees of freedon used in
subroutines DCMPBD and SLVBD (no.
unknowns per node x NN)

node number variable used to build A
vector

soil air porosity

element matric potential I, product of
subroutine DIFFM

air-entry soil matric potential for
element considered, desorption and
adsorption, m H20

average, amplitude of total global

_QRGA

RAV, RAM

RM2

SAA

SAND, SILT, CLAY

SAV, SAM

SIE, SIEA

8M2

SND, SLT, CLY

TAV, TAM

TCON

T81, T31

irradiance, W/m2

matrix used to build EF vector for
moisture

time-varying function for global
radiation

surface sensible heat flow vector
average, amplitude of air relative
humidity

time-varying function for air relative
humidity

matrix used to build EF vector for
moisture

total soil porosity

time-varying function for air wind speed
% sand, silt, clay of B1 horizon
average, amplitude of air wind speed,
m/sec

air-entry matric potential, m H20, of
tilled soil, desorption and adsorption
matrix used to build EF vector for
temperature

% sand, silt, clay of AP horizon
average, amplitude of air temperature,°C
thermal conductivity, product of
subroutine THCON

nodal temperature and moisture content

from previous time step, used as

T32, THZ

TESO

THI, TEI

THO

THSO

TITLE

initial conditions

nodal temperature and moisture content
calculated at present time step
time-varying function for air
temperature

temperature of surface nodes at previous
time step, used for averaging for
surface T's

moisture content and temperature (°C) of
element

moisture content of element's nodes from
two time steps back, used in subroutine
DIFFM to identify whether element is
adsorbing or desorbing

moisture content of surface element from
two time steps back

initial moisture content and temperature
of tilled nodes

vector used to build EF vector for
moisture

moisture content of surface nodes at
previous time step used for averaging
for surface 8's

element volumetric vapor content, 8v,
product of subroutine CAP

immobile moisture content level in soil

title of mainprogram

TO

TOP

TS, THS

ZLEV

ZO

vector used to build EF vector for
temperatures

vector used to build EF vector for
moistures

temperature, moisture content of surface
of surface element considered for
present time step

x and z coordinates of node considered
vector of z coordinates of depths in
grid (surface 2 a 0)

surface roughness, m

ADDITIONAL EXPLANATIONS

Reference

Explanation

 

After statement 28

Comment statement

after statement 35

Comment statement

before statement 750

These are special input statements to
build initial 8, T, and DB gradients

throughout the soil profile

This calculates aVD/az at each
horizontal level in the profile. Soil
profile is first divided into at most
10 horizontal levels. The levels are
numbered 1, 2, 3, etc. from the surface
down, with surface as no. 1. Using air
velocity at bottom level as 0, calculate
the velocity VDU at upper boundary of
the bottom slab. Calculate BVD/az for
that slab and record. Go to the next
slab, above it. Using previously cal-
culated VDU now as VDL, recalculate VDU

for new slab, record new aVb/az, etc.

This calculates BIG/32 for each element
considered. Find the element vertically
in the soil profile between two

recorded levels, assign the avlfiz of

that slab to the element considered.

C T

C EN

1

780

Table B.l. Computer Program
HAINPGH

IME-DEPENDENT BOUNDARY CONDITIONS.
RHAITME.RV.RH)IRVORNMCOSI301516'ITHEIIZ.0-0.ZSII
QRGAITHEoOVuOM)IOV00PMSINI301416*(THE/12o0-0.583))
SAAITHEuSVoSH)SSV*SN‘SIN(3014159.(TPE/IZoO-OoSOBI)
TEAAITHEcTVoTN)3TVOTP‘SINIBolfilbi(TPE/IZoO-Oo583))
DIPENSION 8(3)9C113)cNSI3).EIlZIoQSIIZIoTOI3).TITLEIZO).TDFI3)
DIPENSION THOI3).EEV(3)9AHI3.3IoECH(393)oEHV(3)oECH(3.3)
DIFENSIDN SHZI393IoRP2I3o3DoOHI393)

DIFENSION ZLEV(10)cNLEVI(10)9NLEV2(IDI.DVDDZ(9I

DIPENSION NI(54)9NJ(54).NKI54)gBULK(5#).DVDl(56).AREAISA)
DIPENSION PSIE(SAD-PSIEAISQIoHCONSI5‘).8C(5693o3)oOVSISQD
DIPENSION THSO(60I¢ TESO(40I9THO(40I
DIVENSION SNDISQIoSLTISQIoCLYISA)

DIMENSION XIfiOIoZIQOIoTHIIAOIoTEI(QOIoTHZIAOIoTEZIAOIgAI360)
DIPENSION GGSMI3606I065FI36I06PHII36IOGCHI36936IsGHVI36I
TER GRID DATAoSOIL AND AIR PROPERTIES. INITIAL CONDITIONS.

DATA BULK/5“IOOOeO/QPSIE’5“‘el017/9PSIEA/5‘.’e1017/
DATA IN/l/oIO/BI
DATA SND/Sk'.b35lo$LT/56*o218IoCLY/Sfi‘o1‘7!

DATA THl/40*Oo03/oTEl/‘O‘23o0IoHCONS/Sfi'7ol9E-6/

DATA IDI/O/

READIINoT) TITLE
FORPATIZOAAI
HRITE(ID.780I TITLE
FORMATIIHl/l/IXcZOAA)

C BLILD AN MATRIX.

T96
T95
10

ll

12

13

16

760

00 79‘ I‘Is3

DO 794 J'Ie3

AHIIsJI'OeO

00 T95 I'Io3

AHIIeII3Ie0

READIINOIOI NNsNEsISPsNEVsNNVeNEBsNEHQLEVoITINET.ITQNEVToNIP
FORMATIIZISI

TNIP'NIP

READIIqulI FREQ.AHP.IU,SANC9SILTsCLAYsELEVsTHHP
FORFATIBFIOeSI

JP'NE‘NEVT+A

00 12 J3JP9NE

SNCIJI'SAND

SLTIJI'SILT

CLYIJI‘CLAY

CCNTINUE
PVS‘I1.025..3.SAND‘00026‘.3.SILT‘OecOITT3.CL‘VI‘OeOOI/I10025‘.2.S‘
IND‘O.026MTZTSILT’O.00I“2‘CLAYI

00 I3 N'IVNE

DVSIN)8PVS

READIINpIOI RAV'RAflgTAVgTAHsSAVoSAHsQAVOQAH
FONFATIDFIO-QI

NP'NN

JGF‘NP

JGSPSJGF*Z

JENc=JGSHONP.NDH

DO 760 N‘IoNN

REACIINQITI XINIoZINI

PAINPGN

17 FORMATI2F5e3I
IFIISPeEOeOI GO TO 30
READIINQIBI BLKQSIEQSIEAIHONS'THIQTEI

IO FORMATIbFIOeSI
00 2‘ J'IONEV
RE‘DIINO22I N

22 FCNHATIISI
BULKINI'BLK
DVSINI'OeO252‘EXPI'00192‘TNIPI
PSIEINI‘SIE
PSIEAINI'SIEA
HCCNSINI'HONS

2‘ CCNTINUE
DO 28 J‘IsNNV
READIINsZSI N

25 FORMATIISI
THIINI'THI
TEIINI'TEI

28 CCKTINUE
DU ‘17 J63‘395‘

‘17 BULKIJDI'IDOOeO
DO ‘18 J6=13ela
‘18 BULKIJOI'IZCOeO
DC ‘19 J631902‘
‘I9 BULKIJ6I‘IZSOeO
DO ‘20 J6'25930
‘20 BULKIJ6I3I3‘Oe0
DO ‘21 J6'31036
‘21 BULKIJDI'I‘ASeO
DO ‘22 J6331s‘2
‘22 BULKIJOI'ISTOeO
DO ‘23 J6'9912
TMIIJ6)8.O95
TEAIJOI'ZSeO
‘23 CENTINUE
DO ‘2‘ Jb'laslb
THIIJOI'eIZ
TEIIJOI'ZSeS
‘2‘ CONTINUE
DO ‘25 J63ITQZO
THIIJDI'eI‘
TEAIJbI'ZSeS
‘25 CCNTINUE
DO ‘26 Jb‘IsNN
‘26 THOIJOI‘THIIJOI
DO ‘27 J6'21'2‘
THIIJDI'eIOS
TEIIJOI‘ZbeO
‘27 CONTINUE
DO ‘28 J6325928
THIIJOI'eID
TEIIJDI'ZTQO
‘28 CONTINUE
DO ‘29 J6329s32

FAINPGN

THIIJ6I=e20
TEIIJ6)=ZToO
429 CCATINUE
DC All J6=33o‘0
TEIIJOI=ZT.0
THIIJ6’3020
fill COATINUE
C CONNECT NODES TC ELEMENTS.
30 DO 761 I=loNE
TéI REAO(IN929) NIII).NJII|.NK(II
29 FORHATI3I5I
DO 740 I=IgLEV
740 REAOIINg74ll ZLEVIIIaNLEVIIIIcNLEVZII)
741 FORMATIFIOo‘ozISI
C DLILD BC MATRIX.
DC 788 N=IeNE
NI'NI(N)
NZ‘NJIN)
N3'NKIN)
AREA(NI=(X(NZI‘ZIN3)+X(NII‘ZIN2I+XIN3I¢ZINII-XINZI‘ZINI)-XIN3I’Z(N
IZ)'X(N1).Z(N3II'O.5
BII)‘Z(N2)-Z(N3)
BIZI'ZIN3I-ZINI)
8(3)81(Nl)-Z(N2)
CIII)=X(N3I-X(N2I
CIIZ|3X(Nl)-X(N3I
CII3)=X(NZI-X(NII
DC 50 18193
00 50 J8193
50 BCIN.I0JD8BIII‘B(JI§CIIII‘CIIJ)
788 CCNTINUE
C START DO LOOP FOR TIME STEP.
DO 200 ITINE'IToITIHEToIT
DO 35 I'IvJEND
35 AIII'OoO
TIFE=ITIHE
TIPE‘TIHEI60.0
TIT'IT
TIT‘TIT/bOoO
RHIRHAITIHEpRAV'RAHI
TEAITEAAITIHEoTAVsTAPI
SASSAAITIMEsSAVvSAHI
QRG*0R6A(TIMEoQAVeQAPI
IFICRGeLTeOeOI QRG'OeO
AHP‘AFP‘SA/ISAV+SAHI
C CALCULATE DVDDZ FOR EACH LEVEL.
VCL=0.0
KSLEV-I
DO 745 I’vi
JI’LEV-I
JZ’JI*I
NUPI‘NLEVIIJII
NUPZ'NLEVZIJII
NLOI'NLEVIIJZI

PAINPGH

NLDZ’NLEVZIJZI
TEJI’ITEIINUPIITTEIINUPZII/2.D
TEJ2'ITEIINLDIIOTEIINLCZII/2.D
DELT'TEJZ'TEJI
IFIDELT.LE.0.DT DELT'D.D
DELZ'ZLEVIJZI‘ZLEVIJII
VCL'SQRTI0.0667‘DELT‘DELZ9VDL’.2’
DVCDZIII:IVDU'VDLI/DELZ
IFIDVDDZIII.LT.D.DI CVDDIIII'D.O
VCL'VDU
7‘5 CCNTINUE
DC 811 I'IvNN
GHVIII'D.0
DD DID J‘IoNN
GCPII.JI80.D
DID CCNTINUE
DII CCATINUE
C START DD LOOP FDR EACH ELEMENT FDR TEMPERATURES.
DD 80 N’AQNE
NI'NIINI
NZ'NJINI
N3'AKINI
C CALCULATE DVDDZ FOR THE ELEMENT.
lAVC'IZINII02IN2I’ZIN3II’3.D
K=LEV-I
DD T50 I3I9K
JI'LEV‘I
JZ'JIAI
IFIZAVG.LE.ZLEV(JZIeANDelAVG.GT.2LEVIJIII DDZ'DVDCZIII
750 CCATINUE
TH'ITHIINIIOTHIINZIOTHIIN3II/3.D
TE‘ITEIINII‘TEIINZIOTEIIN3II/3.D
S'I.D'DULKINI/2650.D
P's-TH
C CALCULATE THERHAL CDNDUCTIVITV. HEAT CAPACITY.
CALL THCONITHsSsP9T595NDINI.SLTINIQCLYINI.TCDNI
CALL CAPIS.P.TH.TE.CG.C.THVI
C CALCULATE VELOCITY GRADIENT.
CALL VELIDDlsZINlIsPoDVSINI.FREDOAHPoDVDzIN.)
IFIN.GT.NEDI GO TO ‘8
C CALCULATE EVAPORATION AND HEAT FLOW AT SURFACE.
TS‘ITEIINIIOTEIINZTI/2.D
THS'ITHIINIIOTHIINZII/2.D
THCS'ITHDINIIATHDIN2I§THDIN3II/3.D
CALL DIFFHIS.PSIE(NI.PSIEAINI.THDS.TE.P.SNDINI.SLTINIQCLVINI.TH.TH
IMP.HCDNSINI.HCDNsDIFPsPSII
CALL EVAPIRHOTEA.ZO.ELEV.SA.DRG.PSIOTSQEINI.THS.QSINI.NI
CALCULATE ELEMENT CDND. PATRIX.HEAT VECTOR FDR TEMP.ECUATIDNS.
‘8 PISTCCN‘TIT‘150.D/AREAINI
DZ'C‘AREAINI/Z‘.D
DD 61 I319}
DD 60 J‘Iv3
ECVII.JI33.D.PZ‘DCIN.IoJIOD.D‘QZ‘AHIIoJI
SHZII.JI33.D*PZ‘DCIN.IsJI‘D.D‘DZ.AHII.JI

n

60
61

TO
71

67
72

820
821
C INS

75
76
ac
ea

660

670
675

NAINPGH

CCATINUE

CCATINUE
TOIII‘TEIINI)
TOIZI'TE1INZT
TCI3I=TE1IN3D

DO 71 1=1v3

EF-0.0

DO 70 J'193
ST'SHZII.J)PTC(JI
EF‘EF+ST

CONTINUE

EHVIIifi-EF

CCNTINUE
IFIN.GT.NE8) GO TO 72
DO 67 I8112
EHVII)8EHVIII0(XIN2)-X(N11)‘TIT*1800.0*OS(NI
DO 821 18193
IFIIeEcelI KBNI
IF(I.EO.2) KINZ
IFII.EQ.3I K3N3
GHVIK)=GHV(K)¢EHVIII
DO 820 J8lv3
IFIJ.EO.1) LINI
IFIJ.EQ.2I L3N2
IFIJ.EO.3) L3",
GCPIKoLIIGCNIK9L30ECPIIvJI
CCNTINUE

CCNTINUE

ERT ELEMENT PROPERTIES INTO GLOBAL A VECTOR.
NSIlIlNl

NSIZ)IN2

NSI3I'N3

DO 77 131.3

IISNSII)

JSINP+II
AIJ518AIJ5)OEHVII)
DO 76 J'193
JJ‘NSIJ1-1I+1
IFI44116.16.15_
JS'JGSN*(JJ-1IPNP+II
AIJSIIAIJ5)OECNII.JI
CCATINUE

CCATINUE

CONTINUE
HRITEIIO.84)TIPE
FORMATII/3X.3OHTEHPERATURESICl FOR HOUR NUH. .F6.2//’
DO 660 I'IoNP
GGFIII'AIJGF‘II

DO 675 J'loNBH

DO 670 I'loNP
K8IJ’1)*NP*I
GGSMII.J)8A(JGSH+K)
CCATINUE

CCATINUE

PAINPGH

CALL DCMPBDiGGSF.NP.N8HI
CALL SLVBDIGGSH.GGF.GPHI.NP.N8W.ID)
DD 122 I’loNN
TE2III=GPHIIII
IFITEZIII.LT.TE1III-‘.DI TE2III=TE1III“.O
IFITEZIII.GT.TE1III*‘.DI TE2III3TE1III9‘.D
IFITEZIII.LT.0.DI TE2I1I'D.D
122 CCATINUE
HRITEIIOQZO‘I (I.TE2(II.I'1.NNT
26‘ FORMATIIX.I3.E1‘.5.3X.I3.E1‘.5.3X9I39E1‘.593X.I3.E1‘.5.3X.I3.E1‘.5
1’
DO 120 I‘lsJEND
12D A(I)'0.0
DO 8‘0 I'lgNN
GHVIII30.D
DD 839 J’IONN
CCVIIsJI‘Oeo
839 CCATINUE
8‘0 CONTINUE
C START DO LOOP FOR EACH ELEMENT FOR MOISTURES.
DD 180 N=1vNE
NI'NIINI
NZBNJINI
NB'NKINI
THSITHIINI10THIINZT9TH1IN3II/3.0
TE'ITEIINII‘TEIINZ’PTEIIN3II/3.D
THCZzITHDINII§THOIN2I§THDIN3II/3.D
$31.0‘DULKINI’2650.O
P3S’TH
C CALCULATE VAPOR . VELOCITY GRADIENT. MCISTURE DIFFUSIVITIES.
CALL CAPIS.P.TH.TE.CG.C.THVI
CALL DIFFMISQPSIEINIOPSIEAINI.THDZ.TE.P.5NDINI.SLTINI.CLYIN’.THgTH
1HP.HCCNS(NI.HCCNsDIFPgPSII
CALL DIFFTIHCCN.PST.S.TH.TE.P.DIFTI
C CALCULATE ELEMENT DIFFUS. MATRIX AND EVAP. VECTOR FDR MOISTURE EDUAT.
Sl'TIT“50.D.DIFT/AREAINI
Tl'DIFM’TIT“5D.0/AREAINI
DC 131 I31v3
DO 130 J3193
EDFIIoJI‘TZTDCIN.IOJI‘I‘.0.AREAIN1’12.DI‘AMIIOJI
RM2II9J13T2‘8CIN.I9JI'I‘.D.AREAINI/12.DI*AM(I.J1
QMII.JI=SZ*8CIN9I.JI
130 CCATINUE
131 CCATINUE
TOFIIT'TEIINIITTEZINII
TOFIZI=TE1IN2ITTEZIN2I
TDFI313TE1IN3I§TEZIN3I
THC‘II'THIINIT
THCIZI'THIINZI
THCI31’TH1IN3I
DO 135 13193
EF'Oeo
DO 13‘ J3193
ST'RMZIIpJI‘THCIJI*OPIIoJITTOFIJI

FAINPGH

EF'EFOST
I34 CCNTXNUE
EEV‘I.3’EF*THV.DVDZ‘N"tlt‘120000‘AREA(N'
135 CCKTINUE
lF(NoGToNEB. CC T0 [‘7
DO 138 1'192
138 EEV‘I’3EEV("-(X‘NZ’-X(NI”‘VIT*180000.E(N,
1‘7 00 851 13193
IF11oEQoID KSNI
IFIIoEQoZ) KIN?
IF‘IoEQo3) K=N3
GHV‘K’=GHV(K.9EEV(!’
00 850 J=lo3
tF‘JOEQOl, L‘Nl
IF(JoEQoZ, L=NZ
IF‘JOEQOB, L3N3
GCV‘KoL’36CHLK0L,§EDV[loJ’
850 CONYLNUE
851 CCN‘INUE
{NSERI ELEMENT PROPERTIES INTO A GLOBAL VECTOR.
NS‘I’aNl
’ NS‘Z)‘NZ
NS‘3’3N3
DO 152 1.103
ll‘NSlI’
JS'NP*“
A‘J5"A(J5’*EEV‘I’
DU 151 J8lo3
JJ'NS‘J”IIOL
tFlJJ)15191510150
150 JS'JGSV*(JJ-l"NP*Il
A‘J5’3A(J5.§EDV(IOJ.
151 CCNYINUE
152 CCN‘INUE
180 CCA'INUE
HRIVECIOolefiiTIFE
18‘ FOR'AV‘II3XI2§HHOISYURES FOR HOUR NU". .F6.2/Ii
DC 680 x319NP
680 GGF(["A(JGFOl’
00 685 J319NBH
DO 632 I’lvNP
K3(J'l'.NP‘l
GGSVI[9J’=A‘JGSH§K‘
682 CCK'INUE
685 CCNTINUE
CALL DCMPBD‘GGSNoNPQNBH’
CALL SLVBD!GGSV!GGF.GPHIoNPgNBH.'0’
DO 187 I’loNh
THZ‘I,'GDH[(I'
IFIIHZII).LT.THHPIZ.O) YHZCIOITHHPIZ.O
187 CCAIINUE
HR!‘E(I09265’ (‘0TH2‘L’OI'1'NN’
265 FURNAY(IX9I3951§0593X913151§0593X9l3951‘.5v3xvi3151‘0593‘9‘39E1405
l,

FAINPGH.

DD 190 K‘IONN
TEIIK"TE2(K’
THO‘K’3THICK)
THICKD‘THZCK’
190 CCNTINUE
200 CCN'INUE
STCP
ENC

DlFFH

SUPROUTINE DIFFMIDSqCPSIEoDPSIEAoTHFoDTEoDP'DSAgDSIoDCLoDTHoDTHHPg
lDHCCNS,DHCCN.CCFM.DPSl|
IFIIDTH-THM).GT.0.0) GC T0 £
PSIA=DPSIE
6C T0 6
PSIA=DPSIEA
Dl=CSlIDSA
0235.91*DCL/(DSA+DCL)
0386.2‘SQRTIDID-CZ
BET=Z.619*01**0.2822#(03*0o7l‘*000625#03*‘00125*(CZ+loli“0.0625
DPSCTHs-BfittPSIA‘DS*’EETICTH**(BET+1.0)
H=€PSIIICI.O-BETD)*1DS-DS*‘BET*OTHHP#.(I.O-BET))‘loooo
N8ABS(H)
DHCCN=DHCCNS‘(DTH/DS)“(0.0IS*H+3¢OI
DCFPL‘CHCON*CPSDTH
DCFFV=tlo467E-17‘DP¢((DTEOZTBo0)‘*2o3!‘0PSOTH!/DTE
ODFF‘DDFHL*DCFFV
DPSI8PSIA*(DSIDTH)*’BET
REIURN
ENC

EVAP

SUEROUTINE EVAPIHoTAvDZCcDELEVQDSAOCORGoDPSoOTSoEVPoDTHSoDCS¢NNN)
DATA IN/l/oIC/3/
RI‘9.81’(DELEV-DZO’.(‘A-DISi/ICTA+Z73016).DSAfi‘Z)
[F‘RIOGEOOOI,R180009

513l00/(100'1000.R[’
RA‘(ALOG(DELEVIDZO’)"2/I0¢168.DSA’

RC'RA‘ST

EC'loB‘TA032oO

P5300060“(100365’*.Ec

HA=7o3bth$PSI(TA+273016’
"0‘1.323‘EXP‘11027.0t5/‘23703.015"I'273016.0Ts,
HS'FU'EXPICPS/(46.97.(DTS*273016’Ii

EF'IHS-HA)/RC

DE'2‘5115000*E"

EVP'EF/lOOO-O

0A8125700‘CDT5-‘Ailflc

IF‘O‘OLTOOOO’ QA'OOO

IFCOAoGTo30oO’ OA‘3000
QRL'5067E-8'(TA0273.163“4*(0.60500o048.SORT(1370.0.HA!l
EVIS3O.9O*0318‘DTHS

[F‘CTHSOGTOOOZS’ AL'OOIO

IF‘DTHSOLTOOOIO’ AL'O.25
IF!CTHSoGEoOolO-ANDoCTHSoLEoOoZS, AL=0.l00(0.25-01H53
ORR.(l.0-AL)‘DCRGOQRL-EHIS*5067E-B‘CDTS*273016’#*4
[FCCRNoLToOoO) DE‘OoO
IF‘ORNOGt.OOOO‘NDODEOGVOORN.OO1S’ DE'QRN.OO7S
DQS'QRN-DE‘QA

iF(ORNoGToOoOoAND.D°SoLYoOQOI 005.000

lF!NNN-GT.I) GO '0 902

HRIYE!1099003 ORNoOAoORLvOQSvDE

9CD FORHAT‘IZXo6HCRN ' 9512.392X15HQA ' 9512.392X96HQRL = cElZo3/ZX95H
10$ ‘ .EIZ.3.ZX.6HE ' 961203.
9C2 RETURN
ENC

4C0
410

DIFFT

SLBROUTINE OIFFTCDHCCNcDPSI90$.DTH.CTE.DA.DDT)
DOILa-2o09E-3*DHCCN*CPSI

THLK=00‘O’CS

IF!DTH.GT.THLK) GO TO #00

F=CS

GA'IQO

GO TO #10

F=CA#(1.0+DTH/(DS-THLKl)
GAsoA/(DS-IHLK)
DA1M=5o801E-ll*lDTE+273oOl*'2o3
DTADT=Oo2#*DTH-loSB’CS+2.46
058!lDAODTHtcA)/(0A+O.124’DTHIGA)!*‘2
DDIV=F*DATM*1.024*1.05E-3'0TADT/l.083
DDY=DDTL*DDTV

REYURN

ENC

VEL

SUEROUTINE VEL(DDDZ.ZZ.OA.DCVS.DFR.CAHP.DDVDZD
DK‘S=DDVS**2*l.lllEb‘DAti3l(loO‘DA)“Z
EX=SORTIDFR¢DAI(20.70*DKASD)*12
DVPDZ=CDFR*DAMP/10.35’*EXP(-EX’

480 DDVDZ=DVPDZ+DDDZI§000

‘75 REILRN
ENC

#50

451

THCON

SUEROUTINE THCCN!DTH9CSvDAvCT590519CSIoDCLvDL)

THLK80.§O*DS

SC'AOO.DS

PSAN‘SU‘DSA

PSIL'SO’OSI

PCLAISO‘DCL

02'027Q

52.0611

c2'0521

DLV'I.02§‘1.USE-3*5oBOIE-ll‘iDTEOZ73oOI"2o3*2045166
DLAV'002577E'1*103‘DLV
DKIV80o333‘I2.0]!0.971-0.2§8*0THIOS’0100/(0.22600049bfiDTh/CS”
0L8(DTH‘.595*DKAV‘DA‘DLAVOOZ‘PSAN‘80799OSZ‘PSlL‘2o1790C2‘PCLA*Zo97

15)l(DTHOOKAV‘OA+CZ*PSANOSZ‘PSlL‘CZ‘PCLAD

REYURN
ENC

CAP

SUBROUT‘NE CAPLDSODApoTHoDTEcDCGOCCQDTHV’
DCG330771‘0TE.*19015*10257E3‘(1.0-0.905-6‘DVE*‘10015)
DC'(100’DS"2.011£6*CYH‘§.1956*DA‘DCG
DTVV'O.9OE“6’DTE*’1¢OI5

RE'URN

ENC

DCMPBD

SUEROUTINE DCMPBDIGSFcNPcNBH)
DIPENSION GSP(NP9NBH)
NP1=NP-l
DO 226 lilcNPl
MJ'IONBH-l
IF‘MJOGVONP, MJ‘NP
NJ‘IOI
HKINBE
IF!!NP-I*lioLT.NBH) FKtNP-I+l
N080
DO 225 JSNJ'HJ
MKSPK-l
NDaNcol
NL'NO‘I
DO 225 KslgflK
NK-NDOK
225 GSECJoK)-GSH(JvKl-GSFCIoNLD’GSHIloNKD/GSHIIoI)
226 CChTINUE
RETURN
ENC

SLVBD

SLEROUTINE SLVBD‘GSFoGFoPHIoNPoNBHo[CT
DIPENSIDN GSHTNPoNBHToGFCNpigPHITNPT
DATA TN/l/olO/3/
NP18NP-l
C DECOMPOSITION OF THE COLUFN VECTOR GFI T

DO 250 I‘loNPl
MJ=I+NBH-l
IFTNJaGToNP) PJSNP
NJII+1
L31
DU 250 J8NJ0HJ
L'L‘l

250 GFTJTathJi-GSNTIoLT*GF(IT/GSHTI.1|

C BACKHIRD SUBSTITUTION FOR DETERMINATION OF PHI! T

PHTTNPT=GFTNPTIGSNTN991|
DC 252 K=19NP1
ISNP-K
HJ=NBh
IFTTIONBH-IT.GT.NP) FJ'NP-T*1
SU'3000
DU 251 J=ZoMJ
N=I+J¢1

251 SUFSSUM+GSMTIvJ)*PHI(NT

252 PHITI)8TGFtlT-SUFl/GSFTI'1T
RETURN '
ENE

APPENDIX C
CONVERSION TABLE
The table lists pertinent metric-metric and metric-English

unit conversions.

 

APPENDIX C

Table C.1. Conversion Table

 

 

 

 

 

To get<— Divide by 4 If you have
If you have -——————£>vMultiply by }> To get

bar 0.987 atmosphere
bar 1017 cm water
bar 100 joule/kg
calorie/cm3 . 4.18686 joule/meter:3
calorie/gm 4186 joule/kg

cm water 0.0009703 atmosphere
gram/cm3 1000 kg/meter3
inch mercury 0.03386 bar

joule 0.2389 calorie
langley 1.0 calorie/cm2
langley/minute 697 . 5 watt/meter2
millibar 0.01017 meter water
watt/meter2 2.3898-5 calorie/sec.cm 2

2
watt/meter 86400 megajoule/metergday

 

APPENDIX D

CALCULATION OF METEOROLOGICAL INPUTS

Explained here is the process by which the diurnal cycle
inputs for air relative humidity, air temperature, air wind

speed, and total global short-wave irradiance were calculated.

APPENDIX D
Calculation of meteorological inputs.
Air temperature, global irradiance and wind speed all were
,modeled in the following form:
v - M + A Sinv n (t/lZ - k)
and relative humidity as the form:
V - M + A Cos n (t/12 - k)

where V - meteorological variable

daily mean
amplitude

‘time (hrs)

3‘")!
a

- phase shift to peak in early afternoon

Air temperature, wind speed and relative humidity data
were found in the literature in enough detail where M and A
could be calculated directly. However, global irradiance is
normally published as daily, or even monthly, averages.

To fit a diurnal cycle to a daily average figure for
radiation (see Figure D1), area A was set equal to area 8
using integral calculus. Published values by Bruce, et a1.
(1977) and Van Bevel and Hillel (1976) indicated that area 8
was larger than that under half a sine wave, thus QAV was
set equal to a value greater than 0, and the sine wave for

radiation was truncated below the 0 line. Thus:

18 set
I [QAV + QAM Sin w (t/12 - 0.5)] dt =

6
. Q(published daily avg.) x 24 hours

The integration produces one equation with 2 unknowns QAM

and QAV:

l2 QAV + 7.64 QAM - 0 (published) x 24
Fitting to published data, the author chose QAV - 50 W/m2.

and used the equation to calculate the corresponding QAM.

published 0 avg.

P/ 4/ 9| 53

F 24 hrsf)‘

 

 

 

 

 

 

Figure 0.1. Global irradiance cycle computation.