Modeling coupled water-heat flow and impacts upon chemical transport in mulched soil beds


Material Information

Modeling coupled water-heat flow and impacts upon chemical transport in mulched soil beds
Physical Description:
xii, 153 leaves : ill. ; 29 cm.
Shinde, Dilip, 1962-
Publication Date:


Subjects / Keywords:
Soil management   ( lcsh )
Mulching   ( lcsh )
Soil chemistry   ( lcsh )
Soil and Water Science thesis, Ph.D   ( lcsh )
Dissertations, Academic -- Soil and Water Science -- UF   ( lcsh )
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph.D.)--University of Florida, 1997.
Includes bibliographical references (leaves 142-152).
Statement of Responsibility:
by Dilip Shinde.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 028638406
oclc - 38863565
System ID:

This item is only available as the following downloads:

Full Text






I dedicate this work to my parents, father (late) Krishna Rao Shinde and

mother Vatsala Shinde, who instilled in me the importance of education and



I express my sincere gratitude to Professors R.S. Mansell (chair) and A.G.

Hornsby (cochair) for their guidance and support during this research work.

Appreciation is also extended to the members of my advisory committee: Professors

P.S.C. Rao, B.L. McNeal, and F. S. Zazueta; for useful comments, help and support.

I would like to extend my appreciation and thanks to Dr. B.P. Mohanty (Univ.

of California, Riverside) for providing numerous technical literature citations.

Thanks are also due to Professor P.N. Kizza for encouragement and support. I would

like to also thank Dr. A. Akpoji and Ron Jessup for useful discussions on finite

element methods and numerical programming, respectively. I acknowledge Dr. John

L. Neiber (Agric. Eng. Dep., Univ. of Minnesota) for providing the FORTRAN code

used extensively in creating finite element-mesh in this work.

I take this opportunity to thank also 'The Rotary Foundation of The Rotary

International' for the prestigious "FREEDOM FROM HUNGER" scholarship to

enhance my knowledge and join the fight to free the world from hunger.

I would like to acknowledge the sit i- r a r, i r1, r ,:.- l ':. I .rT provided by

Professor A.G. Homsby during the latter part of this research. I would also like to

acknowledge the facilities and partial financial support provided by Professors R S.

Mansell and B.L. McNeal, respectively, during this research.

Lastly, appreciation is extended to my family for their help and support.



ACKNOWLEDGMENTS ........................ ................ iii

LIST OF TABLES .................... ............ .............. vii

LIST OF FIGURES ......................................... .......... ix

ABSTRACT ............ ................................... xi

INTRODUCTION .................. ........................ ......... 1

Coupled Versus Uncoupled Flow Processes in Soil ..................... 1
Soil Therni F-ri .......................... .... .. .2
Non-isothermal Soil Conditions ............................ 3
Mechanistic Mathematical Modeling .............. .... ...... 4
H ypothesis ..................................... ........ 6. 6
Objectives ....................................... ......... 7

REVIEW OF LITERATURE .................................... .. ..... 8

Isothermal Water Flow in Porous Media ............ ............ 9
Soil Temperature and Heat Flow ............... ............. .11
Soil Thermal and Water Regimes under Mulching ................... .. 15
Coupled Heat and Water Flows ................................... 19
Philip-de Vries Theory (PDT) ................................ 20
Application of Philip-de Vries Theory .............. ... ...... 22
Solute Transport during Coupled Heat and Water Flow ................. 24
Temperature Effects on Solute Transport ........................... 27
Summary ...............................................29

MATHEMATICAL MODEL ............................................. 31

Liquid and Vapor Water Flow in Soil ................... .. ..... 32
Heat Flow in Soil ................ .. .............. .......... 32
Transport of Multiple Species of Solutes in Soil ....................... 33
Energy Balance at the Soil/Atmosphere Interface ................. 34

Initial and Boundary Conditions .................................... 38
Input Parameter Requirements ................................. 39

MODEL VERIFICATION AND VALIDATION ............................... 44

Isothermal Water Flow ............... .................. ...... 45
Uncoupled Heat Flow ................ ......................... 46
Coupled Water and Heat Flow ........................ ........ 47
Isothermal 1-D Solute Transport during Steady Water Flow ............. 48
Isothermal 2-D Solute Transport during Steady Water Flow ............. 49
Isothermal 1-D Solute Transport with Nitrification Chain during Steady
W ater Flow ............................ ...... ......... 51
Isothermal 1-D Solute Transport involving Gas Flow ................... 52
Summary ............... ....................................54

MODEL APPLICATION ................ ............................... 55

Fate and Transport of Methyl Bromide Fumigation in Soils ............. 55
Transport Mechanisms for MBr Fumigant in Soil ................ 59
Fate of MBr Fumigant in Soil .............................. 60
Permeation/Transmission of MBr through Polyethylene Mulch .... 63
Modeling o! IE. TI in ri.:.. n .iuriri. Fumigation ................ 68
Analysis of Fate and Transport of MBr Fumigant .................... 72
Approach ........ ................................... 73
Methodology ................... ............ .. ...... 75
Results and Discussion .......................................... 79
Plastic-Mulching of Soil Beds to Control MBr Loss to the
Atmosphere ............... .......... ... ...... ... 79
Non-Isothermal versus Isothermal Soil-Conditions .............. 85
Variable Water Saturation of Soil .......................... .. 87
Depth of Fumigation Injection ................... ............ 97
Effect of Irrigation (Water Dousing of the Bed Surface) ........... 98
Fumigation Efficiency Analysis ............................. 101
Hydrolysis of MBr (Sensitivity Analysis) ....................... 107
Diurnal Water and Heat Dynamics in a Plastic Covered Soil Bed .. 108
Summary and Conclusions .................................... 119

APPENDIX A SOIL MOISTURE FLOW .................. ............ .. 122

APPENDIX B SOIL HEAT FLOW ................. ........... .. 125

APPENDIX C CHEMICAL TRANSPORT .................. ........... .127

APPENDIX D NUMERICAL APPROACH ................................ 132

LIST OF REFERENCES ................... ..................... .... 142

BIOGRAPHICAL SKETCH ............................................. 153


Table page

2-1. Average thermal properties of major soil constituents at 20 OC and
1 atm (Ghildyal and Tripathi, 1987) .............................. 12

2-2. Short-wave optical properties of five plastic mulches used in commercial
crop production. Data adapted from Ham et al. (1993) .................. 17

5-1. Earlier estimates for source-contributions of MBr to the atmosphere
(Grojesan, 1991; Khalil et al., 1993; Lobert et al., 1990; Mand and
Andreae, 1994; Yagi et al., 1993; and Zurer, 1993) ................. .... 56

5-2. A recently revised budget for atmospheric MBr (Yvon-Lewis and
Butler, 1997) ......... .. ................. ... ... ...... 57

5-3. Selected chemical and physical properties of MBr (CH3 Br) (Adapted
from Hacherl, 1994) ................... ..................... 58

5-4. Measured degradation rate coefficient and half life of MBr as a function
of soil water content (Adapted from Jin and Jury, 1995) ................ 63

5-5. Rate of escape of MBr through several films at various temperatures
(Adapted from Kolbezen and Abu-El-Haj, 1977) ................... .. 65

5-6. Phi,: i,_J Properties for 5 distinct honzons of a profile of Arredondo fine
sand soil .... ................... ..... ................... 77

5-7. Parameter values for 5 distinct horizons of a profile of Arredondo fine
sand soil ................ ....... .............. ... ........ .. 78

5-8. Material balance (unit: %) of MBr under measured and enhanced levels
of diffusion coefficient for plastic mulch ........................... 84

5-9. Material balance (unit: %) of MBr for fumigation under isothermal and
non-isothermal conditions ..................................... .86

5-10. Material balance of MBr (unit: %) under different water saturations for
one and 3 days after furm gant injection ........................... .. 90

5-11. Matenal balance of MBr (unit: %) during fumigation under different
water saturations for 65-cm injection zone .......... .......... 98

5-12. Material balance (unit: %) of MBr under water dousing and no-dousing 100


Ficrure pEae

2-1. Measurements of soil temperature at 10-cm depth beneath five different
plastic mulches and a bare plot near Manhattan, KS on day 195 of the
year 1991. Adapted from Ham et al. (1993) ................... ....... 18

4-1. Observed and computed distributions of water content in a sand column
during constant flux infiltration ................... .. ......... 45

4-2. Observed and computed temperatures for two depths in a field soil ........ 46

4-3. Observed and computed water content after 31 days under the influence
of both matric head and thermal gradients .................. .. 48

4-4. Observed and computed breakthrough curves for boron ................ 49

4-5. Location of solute concentration fronts for three times as estimated
numerically and with the analytical solution .................. ...... 50

4-6. Analytically determined and numerically computed concentration profiles
at 200 hours ................................... ............... 52

4-7. Analytically determined and numerically computed concentration profiles
at 1, 3, and 5 hours ............ ......................... ...... 53

5-1. Temperature dependence of MBr diffusion in water (Data from Maharajh
and Walkley, 1973) ........................................... 60

5-2. Laboratory-measured and estimated Henry's Law constants at different
temperatures for partitioning of MBr into water (Adapted from Mutziger
et al., 1996) ................... .................. ...... ..... 61

5-3. MBr emission fluxes in plastic mulched beds during fumigation at the
University of Florida Green Acres Experimental Station located near
Gainesville, FL ................ ............................ 66

5-4. MBr diffusion through black LDPE of different thicknesses (Data from
Hacherl, 1994; and Kolbezan and Abu-El-Haj, 1977) ....... ........ 68


5-5. Simulated and measured flux density through plastic-mulched surface.
Adapted from Wang et al. (1997) .................. ............ 71

5-6. Schematic of the discretized simulation domain of the soil bed ........... 76

5-7. Measured and elevated levels of diffusion coefficients through 1.25-mil
black LDPE plastic-mulch used in model simulations .................. 80

5-8. Losses of MBr to the atmosphere through soil bed and furrow surfaces
using measured and elevated diffusion coefficients for MBr through
plastic mulch ...................................... ......... 81

5-9. MBr total concentration (pmol cm3 ) contours and gaseous MBr flux
vectors (pmol cm2 hr ') in the bed 6 hrs after injection of MBr fumigant... 82

5-10. Emission fluxes of MBr through bed and furrow under isothermal and non-
isothermal transport of fumigant ............................... 85

5-11. Emission fluxes of MBr through bed and furrow during fumigation
performed under different profile water-saturation scenarios ............ 88

5-12. Contours of water content (cm3 cm-3), temperature (C), and MBr
concentration (umol) with gaseous fluxes (pmol cm"2 hr ') of MBr under
2 different soil-water saturation scenarios after fumigation ............. 91

5-13. Emission fluxes of MBr under different water saturations for a 66-cm
injection zone during the first 3 days after fumigant injection ........... 97

5-14. Effect of water dousing on the MBr emission fluxes from the soil bed
during fumigation ................................................ 99

5-15. Contours of water content, temperature, and MBr concentration with
gaseous fluxes of MBr at 9 hours after fumigation under water dousing
and no-dousing on the bed-surface ....... ........................ 100

5-16. Soil sterilization zones for different scenarios 3 days after MBr fumigation
using different criteria ......................................... 104

5-17. Loss of MBr (unit: %) under different rates of hydrolysis of MBr ......... 107

5-18. Diurnal water and heat dynamics in a mulched bed for a wet soil. ....... 113

5-19. Diurnal water and heat dynamics in a mulched bed for a dry soil. ....... 116

D4-1 Linear triangular element ........... .............. .. .. 136


Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor o! F r~l..:;.:i ri.



Fl l .- i L r- .L,

December 1997

Chair: Robert S. Mansell
Cochair: Arthur G. Homsby
Major Department: Soil and Water Science

Plastic-covered soil bed systems are widely used in the southeastern United

States for vegetable and fruit production. Plastic-mulching results in highly variable

and transient soil thermal regimes both in space and time under atmospheric

influences (e.g., solar radiation, wind, etc.). Much evidence exists that water flow is

affected by non-isothermal conditions and vice-versa. In this work, a two-

dimensional model for simultaneously describing chemical transport and coupled

water-heat flows in plastic-covered soil beds was developed with many

improvements over earlier efforts. Realistic plastic-mulch boundary conditions are

thoroughly described by including optical properties of the plastic mulch into the

energy balance at the soil/atmospheric interface along with other atmospheric

influences of weather components. This model provides opportunity to consider the

effect of coupled water-heat flows upon estimates of chemical transport (both in the

liquid and gaseous phases) in plastic-mulch culture production systems.

Chemical transport in soil beds in the model utilize consideration of

chemical existence in solid, liquid, and gas phases. Temperature dependence of the

rates of chemical reactions and volatilization is described efficiently by

incorporating temperature dependence in fate and transport coefficients. The model

uniquely provides opportunity to simulate coupled water-heat flow and multiple-

solute transport with a complex plastic-mulch atmospheric boundary for realistic

field situations.

The model was used to analyze the fate and transport of methyl bromide

(MBr) during fumigation of mulched soil beds. MBr, a highly volatile pesticide, is

widely used in California, Florida, and Hawaii for its broad-spectrum efficiency in

controlling insects, fungi, and weeds. Unfortunately, MBr is a major environmental

hazard due to its depletion effect on the atmosphere ozone layer. Model analyses

revealed that non-isothermal conditions should be considered in the accurate

description of MBr fate and transport. Modeling results indicate that utilizing a

relatively wet soil profile, deeper injection, and water dousing on the bed-surface

before installation of plastic mulch will decrease atmospheric emissions of MBr.

However, deeper injection of MBr fumigant, although environmentally less

hazardous, may decrease effective sterilization of crop root zones in the soil.


Coupled Versus Uncoupled Flow Processes in Soil

Flow of water, air, and heat in soil is generally assumed to be uncoupled,

that is, the driving forces causing individual flow are generally considered to be

mutually independent. Darcy's equation for liquid water flow, based upon spatial

gradients of soil-water potential, and Fourier's equation for heat flow based upon

spatial gradients of soil temperature, provide two examples of uncoupled flow

processes. In nature, however, water is also known to move in response to

molecular forces, such as thermal and osmotic gradients. Similarly, considerable

heat flow may occur as latent heat by water vapor diffusion in the air-filled soil

pores, and as sensible heat transfer due to liquid and vapor movement. In a partially

saturated soil, flow of water due to thermal gradients occurs through zones

occupied by liquid as well as air, and is assisted by multiple condensation and

evaporation cycles Pri, p aip de Vries, 1957). Thus, transport processes occumng

in the soil under field conditions are intricate and complex.

Diffusion theory indicates that water molecules are in constant thermal

motion, accompanied by jump-like transitions from one temporary equilibrium state

to another (Ghildyal and Tripathi, 1987). The frequency of transition and the number


of molecules traversing a certain section of the liquid is highly dependent upon the

temperature gradient and causes self thermo-diffusion of water. As temperature

gradients are accompanied by gradients of surface tension at the air-water

interface, there is a possibility of thermo-capillary flow and thermo-capillary film

flow of water in the liquid phase (Ghildyal and Tripathi, 1987). Thermo-capillary

flow occurs in those capillaries where liquid-water menisci exist. Here, temperature

variation induces the variation of capillary pressure and a corresponding shift of the

liquid column. Thermo-: j 1 ll r, I n i ... is caused due to gradients of surface

tension (change in water affinity with change in temperature) at the water-film

surface. Both of these mechanisms cn.rarrt-uir -i. _nri, arri ,'to the total thermal

water transfer under moderate to high water content conditions. Movement of water

in the liquid state due to changes in the volume of entrapped air with temperature

has also been observed (Ghildyal and Tripathi, 1987). In water-saturated or nearly

saturated soil, the entrapped air may cause direct transport of water through the

pressure it exerts and, in unsaturated soil, through an apparent increase in the

temperature coefficient of the matric potential. Convective flow of water vapor may

also occur due to differences in the densities of cold and warm air.

Soil Thermal Regime

The thermal regime of a soil depends upon radiative solar heat flux into the

soil and upon heat-transfer processes occurring in the soil and between the soil and

air which, in turn, depend on the thermal charactenstics of soil solids, gases, and

water. The soil thermal regime, which greatly modifies the overlying microclimate,


is generally characterized in terms of soil temperature. Soil temperature fluctuations

are especially important since they impact physical, chemical, and biological

processes in the soil.

The source of heat for soil is incoming net solar radiation. On a clear day,

Geiger (1957) described the typical distribution of incoming solar radiation as: 1)

global radiation, which is the sum of the direct-beam penetration to the earth (19 %)

and the diffuse-sky radiation (26 %); 2) scattered radiation returned to outer space

(11 %); 3) radiation reflected by clouds (28 %); and 4) radiation absorbed by the

atmosphere (15 %).

Non-isothermal Soil Conditions

Much published literature confirms that soil-water movement is affected by

spatial gradients of both temperature and water potential. Thus, under non-

isothermal conditions, movement of water (both liquid water and vapor tend to

move from regions of high to low temperature) through soil should be considered for

both liquid and vapor phases. Non-isothermal conditions have broader implications

in that solute-transport mechanisms are associated with water movement, and the

fate and transport of volatile organic chemicals are also influenced by the soil

thermal environment.

The influence of thermal gradients upon water movement in soils has largely

been observed in two areas (Jury, 1973): 1) in the hot, dry surface layer of a bare soil

where the interface between the soil and atmosphere causes water movement to be

largely in the vapor phase; and 2) deeper in the soil where thermal gradients persist

over seasons and can slowly transport liquid upward or downward. In the special

case of plastic-mulched soil, imposed thermal gradients tend to be greater at the

soil surface and the mulch also provides a partial barrier to vapor flow between the

soil and atmosphere and causes a return flow of water after condensation under the


Both water movement and contaminant transport are influenced by

temperature gradients in soils. Natural temperature gradients are usually highest

during daytime (midday) near the land surface, but the relatively small geothermal

gradient below the level of seasonal temperature variation can be a significant

factor in determining the percolation of water in deep unsaturated zones (Evans and

Nicholson, 1987).

Mechanistic Mathematical Modeling

Mechanistic mathematical models provide powerful tools to describe and

interpret flow and transport processes in soil-water systems. Bear and Verruijt (1987)

defined a model as a 'simplified version of the real system that approximately

simulates the excitation-response relations of the latter'. The first step in modeling

involves the construction of a conceptual model of the problem, consisting of a set

of assumptions that reduce the real problem and the real domain -. -- ;r.- i-,i, -e

versions (ultimately permitting mathematical description of the processes) that are

acceptable in view of the modeling objectives. In the second step, the conceptual

model is expressed in the form of a mathematical model (Bear and Verruijt, 1987).

Properly posed mathematical models for initial-value problems involving mass and

energy transport in porous media require (i) appropriate differential transport

equations, (ii) initial conditions in the flow domain, and (iii) boundary conditions

imposed along the surfaces of the flow domain. Models are essential in performing

complex analyses and in making informed predictions. Friedman et al. (1984) stated

that, in some cases, models have increased the accuracy of estimates of future

events to a level far beyond 'best judgement' decisions. Anderson and Woessner

(1992) emphasized that models provide the best way to make an informed analysis

or prediction concerning consequences for a proposed action. The authors further

stated that models can also be used in an interpretive sense to gain insight into the

controlling parameters in a site-specific setting or as a framework for assembling

and organizing field data and formulating ideas about system dynamics.

Schwarzenbach et al. (1993) defined a model in the following way: 'A model is an

imitation of reality which stresses those aspects that are assumed to be important

and omits all properties considered to be nonessential'. These authors summarized

the role of modeling explicitly in the following ways:

1) Mathematical models allow knowledge acquired for different systems and/ or

different situations to be combined in order to provide the ultimate aim of

constructing general theories. Models permit making predictions, which then must

be verified by measurement.

2) Classic scientific theories begin as models, and it is only long-time positive

experience with some of these models which leads scientists to believe that they are

more true and fundamental than others and should thus bear the name 'law'

3) In environmental sciences, models serve to sort out alternative explanations for

observations made in nature which cannot be controlled in the same way as

conditions for a laboratory experiment can be controlled.

4) Models should always be constructed to transfer the knowledge acquired in one

system to another situation, thus enabling a deeper understanding of the processes

responsible for the measured data. A model should be applied in order to design

new experiments or observation programs that are critical for the testing of


5) Real observations are always superior to outputs produced by a given model. But

model results can help us examine various scenarios and likely future outcomes

regarding the behavior of environmental systems for which real data would never be

collected as frequently and as ubiquitously as needed.


The soil thermal regime within plastic-mulched soil beds is well-known to be

influenced by interactions between the overlying atmosphere and the plastic-

covered soil surface. The resulting non-isothermal conditions affect water flow and

redistribution. Hence, coupling of water and heat flows is necessary to model

effectively both components within the vadose zone, especially the crop root zone.

The hypothesis here is that model description of chemical transport within plastic-

mulched soil beds under non-isothermal field conditions requires the utilization of

coupled water and heat flows in order to provide reahstic results for management

practices such as fumigation with methyl bromide to control pathogens. Impacts of

this coupled process upon chemical fate (transformations, microbial degradation,

etc.) and transport are expected to be particularly significant for volatile chemicals

such as methyl bromide that may appear in aqueous, gaseous, and solid phases of

soil water systems. Information obtained here regarding the dynamics of the fate

and transport of methyl bromide during fumigation of agricultural soils should also

be useful in evaluating future fumigation with fumigants other than methyl

bromide, which will no longer be used in this country after the year 2001.


This research was performed to model and investigate the interactive flows

of water, heat, and solutes in non-isothermal soil environments (the vadose zone) in

order to minimize empiricism from existing knowledge of such systems. The

following objectives were formulated:

1) Develop a mathematical model describing the coupled processes of water,

heat and chemical transport (2-D) for mulched soil-bed systems, where solar energy

imposes non-isothermal conditions.

2) Verify, validate, evaluate, and refine the model as needed for cases of both

coupled and uncoupled flows of water, heat and chemicals with existing

appropriate models (analytical) and published data.

3) Use model simulations to determine impacts of coupled water and heat

transport due to non-isothermal conditions upon the transport of a volatile pesticide

such as methyl bromide during fumigation of plastic-covered soil beds


The rate of liquid water movement in soil is important to the study of solute

transport, nutrient availability, efficacy of pesticides, and also water availability to

crop root systems during agricultural production. Physical properties of the soil

influence water movement in two important ways: the size distribution and

structural arrangement of soil particles determine the pore space configuration

through which water moves; and the interaction between the soil and water gives

rise to water-moving forces. In water-saturated soils such as occur below

groundwater tables where pores are completely filled with water, the fluid is a

single phase and forces originating from gravity, including the pressure gradient,

determine the direction and magnitude of flow. The energy potential for the pore

water under saturated conditions is generally positive, or greater than the

atmospheric pressure. The saturated flow process for the simple isothermal case is

described mathematically by Darcy's equation and Laplace's equation.

Generally, under field conditions, soil pores are incompletely or partially

saturated and the energy potential of water is negative (i.e. a suction), being less

than atmospheric pressure. It depends not only upon position in the gravitational

field, but also upon sorptive forces associated with interfacial boundaries. If the

unfilled pore space is filled with gases such as air, movement in the resulting

composite system of solid, liquid and gaseous phases is referred to as unsaturated

flow. This flow condition in soil is particularly important when the groundwater

table occurs at considerable depth below the soil surface (below the vadose zone).

Isothermal Water Flow in Porous Media

Under saturated conditions where the pores of the soil system are completely

filled with water, the transmission coefficient (i.e., hydraulic conductivity) depends

principally upon the viscosity (u,(T)) and density (p,(T)) (both being functions of

temperature, T) of the water and the nature of the porous media. The volumetric flux

of water flow through soil or other porous material is usually given by Darcy's

equation as

q = -K VP [2.1]

where q, is liquid flux, K is hydraulic conductivity, and P is hydraulic pressure of the


Except for the isothermal saturated condition, Darcy's equation, as originally

postulated, does not apply precisely to water flow through soil. However, this

equation has been used since 1856 for describing liquid flow through porous

material with little reservation.

Under unsaturated conditions when the soil contains both liquid and

gaseous phases, the most frequently used modification of Darcy's equation involves

substituting the matric-potential gradient for the hydraulic pressure gradient as the

driving force. Under these conditions, hydraulic conductivity is no longer a constant

but is a function of capillary potential as well as temperature, which yields

q = -K(h, T)VH ; H = h(L, T)+z K= k ) g [2.2

Here His total potential, h is matric potential, z, is gravitational potential (i.e.,

elevation), 8L is liquid volumetric water content, T is temperature, k, is intrinsic

permeability and a. i, tih acceleration of gravity. Under isothermal conditions the

temperature dependence is eliminated and total potential gradient is the driving

force. Richards (1931) derived the equation of isothermal flow through a porous

medium by using the continuity equation .-1..:n ri-ri ni i-. conservation) and

Darcy's equation. Assuming that the medium is incompressible, the continuity

equation for water flow through a homogenous medium was given as

a(p = -VpqL [2.31

where t is time and PL is the liquid density. In very dry soils, the capillary or matnc

potential attains great negative values and, in the unsaturated region above a water

table, the capillary driving force dominates the gravitational force.

Isothermal theory has limited application for soils under actual field

conditions, where heat flow is transient in time and varies in space. It has been

observed that heat is evolved at the wetting front when water infiltrates a very dry

soil of fine texture (Anderson and Linville, 1960). Cooling of soil occurs when water

is evaporated from soil (Wiegand and Taylor, 1962) and warming occurs due to

release of thermal energy during condensation of water vapor. Absorption of solar

radiation at the soil surface and chemical or biological activity within the soil may


also result in temperature gradients. Thus, temperature gradients must be included

in the flow equation in order to effectively describe water flow in surface soil

characterized by transient thermal regimes.

Soil Temperature and Heat Flow

Under field conditions, the profiles of soil temperature change rapidly during

a normal day. When incoming solar radiation falls on the soil surface the soil surface

becomes heated, causing a difference of temperature between the surface and the

subsoil. This difference causes heat to flow downwards as a temperature wave. As

the heat wave penetrates slowly downward in response to soil thermal properties,

moisture content, and structure of the soil, a phase shift occurs in the temperature

curve with depth which leads to a time lag in the occurrence of temperature

maxima at different depths in the soil. Also, the amplitude of the temperature wave

decreases with soil depth due to storage of thermal energy. During night time, when

the soil surface cools down, the lower portion of the soil profile remains warmer

than the air temperature. Annual temperature fluctuations extend much deeper into

the soil compared to diumal changes, and are also in general periodic in nature.

Soil texture has a distinct influence on soil heat flux and, hence, on soil

temperature. The amplitude of the daily temperature wave decreases in the order

sand > loam > peat > clay. This order can be attributed to the greater heat

capacity of the finer-textured materials and their normally greater thermal


Table 2-1. Average thermal properties of major soil constituents at 20 C and
1 atm (Ghildyal and Tripathi, 1987).

Volumetric Heat Thermal
Material Specific Heat Thermal Conductivity Thermal
(cal -10o-) Capacity 3 1 1 I Difusivity
(calg oC) (calm- C-1) (10 cal cm- s- C-) Diffusivit

Quartz 0.175 0.46 20.00 43.00
Other Minerals 0.175 0.46 7.00 15.00
Organic Matter 0.46 0.60 0.60 1.00
Water 1.00 1.00 1.42 1.42
Air 0.24 0.00029 0.062 0.21

Sandy soils i-ri-r.ajl ,' hold less water and usually drain faster than finer-

textured soils. At low water contents, air occupies much of the pore space and

hydraulic conductivity of the soil becomes limited. As such, thermal capacity and

conductivity of soil with a low degree of water saturation is governed by the

volumetric content of water, which has the highest heat capacity (Table 2-1) among

the substances commonly found in soils. Because of very small thermal diffusivity,

air provides thermal insulation against rapid changes in temperature.

Three mechanisms by which heat flow occurs in soils include conduction,

convection and radiation (Hillel, 1980). Radiative heat transfer is of importance only

in dry soils at high temperatures and within large pores. It is the least-important

mechanism, and is often neglected. Convection can be separated into two

components, forced and free. Forced convection is mostly negligible, with the

exception of infiltration during irrigation or heavy rain. Free convection occurs as a

result of fluid density gradients. Convective heat transport in the vapor phase can be

Spair, r:.r-r,. ,, r.:, ein i ,l, and latent heat components. The sensible heat capacity of

soil air is very low (Table 2-1) and, as soil air flow is most often small, latent heat

transport is of primary importance. Latent heat transport occurs when a fluid such

as water changes phase and is transported as a new phase (such as vapor).

Vaporization of liquid water requires heat (540 cal g-') that is transported in the

vapor phase through a vapor-phase gradient and may be released during

condensation at a location different from where the liquid originated.

A very simple theory of heat transfer in soils is based upon uncoupled mass

and heat flows and considers only heat conduction, without regard for any other

transport mechanisms which may actually occur. The thermal flux in the soil

(assuming that thermal gradients do not result in water flow) is defined from

Fourier's first law of heat conduction as

q = -XVT [2.4)

where qh is thermal flux and A is thermal conductivity. Invoking the principle of

energy conservation in the form of the thermal continuity equation (Hillel, 1980)

8T x
S= VDH VT; with D, (2.51

where DH is thermal diffusivity and C, is volumetric heat capacity of the soil.

Conduction of heat in soils occurs through the solid, liquid, and gaseous

phases. Since these phases have different heat capacities and thermal

conductivities, the main complication when solving the heat conduction equation in

soils is caused by determining the average soil thermal properties, de Vries (1963)

reported a method to calculate heat capacity of soil material by addition of the heat


capacities of the various constituents, weighted according to their volume fractions.

The heat capacity as such can be found by adding the heat capacities of the various

phases present as

C = (xsiC s)+X C +XaC [2.6]

where the subscripts s, w and a refer to solid, liquid water, and air phases

respectively; i distinguishes the various materials in the solid phase; and X denotes

the volumetric fraction of a given phase (de Vries, 1963). Thus, volumetric heat

capacity for soil clearly depends upon water content.

The thermal conductivity of a soil also depends in a rather complicated way

on the composition of the soil, in particular its water content. A physical model,

assuming sohd and gas particles to be suspended in the liquid phase, for estimating

the thermal conductivity in soils was given by de Vries (1963) as

x w+ E (k *X)1t +k.*X/
X = X k X +k-X [2.7]

The summation for this equation extends over the different solid soil constituents,

which are characterized by their thermal conductivities and shapes. The

multiplication factor k' in Eq. [2.7] represents the ratio of the spatial average of the

temperature gradient in the soil grains of kind i and the spatial average of the

temperature gradient in water.

The influence of latent heat transfer in air-filled pores is proportional to the

temperature gradient in these pores and can be taken into account by adding the

thermal conductivity of vapor (Lp) to the thermal conductivity of dry air (X,). The

1,, term represents the increase in gaseous-phase thermal conductivity due to

evaporation and condensation of water vapor. At complete dryness and at very low

water contents, the model assumes air rather than liquid to be the continuous

medium. However, at this stage (i.e. completely dry soil), de Vries (1963) proposed

introduction of a multiplication factor of 1.25 to Eq. [2.7] (i.e. 1.25 1). de Vries

reported that errors in calculating X in this manner would not exceed 10 %, unless

the soil grains have a different shape than those assumed.

Wescott and Wierenga (1972) demonstrated the need for coupling heat and

water movement when predicting heat transfer by conduction and vapor movement

for a dry soil. They found that vapor movement accounted for 40-60 % of the heat

flow in the top 2 cm of the soil pr..,' ii- d u rig certain periods of the day and for

20-25 % of the total heat flow even at the 25-cm depth.

Soil Thermal and Water Recimes under Mulchina

Mulches are used in agriculture not only because they suppress weeds and

reduce water loss from the soil surface, but because they also modify the thermal

microclimate of the soil. Rosenberg et al. (1983) presented a good review of mulch

effects on soil processes and defined mulching as the application or creation of any

soil cover that constitutes a barrier to the transfer of either heat or vapor, or both.

Waggoner et al. (1960) applied black plastic, translucent plastic, aluminum

foil, paper, and hay mulches to soil surfaces and measured the resulting alteration

in surface energy balance and in soil and air temperature profiles. Black plastic

reduced the outgoing radiation as was evidenced by the greater net radiation,

whereas hay and paper (and especially paper) increased it. All mulches decreased

the quantity of energy consumed in evaporation by blocking the transport of vapor

out of the soil. Sensible heat generated at the surface was increased by the black

plastic and the hay, and the temperature of the hay and black plastic was elevated.

Black plastic increased in soil temperature, whereas hay and paper had a cooling

effect and decreased the amount of heat penetrating the soil. Ekem (1967) in Hawaii

found that black plastic mulch reduced water use by pineapples grown in

lysimeters. Note here that the well-known practice of surface mulching (for any

purpose) decreases evaporative water loss by creating a layered condition in the soil

due to differential water contents.

The diurnal supply of solar radiant energy at the soil surface during daylight

hours and loss during night-time hours is the major determinant of the thermal

microclimate of soil. The presence of mulch material at the soil surface creates a

-1,. rr,i ,.nfr. effect (similar to the effect of soil layers within the soil). The use of

nonporous mulch such as polyethylene plastic also creates a barrier to water vapor

flow to the atmosphere.

Optical properties of the plastic mulches play an important role in interaction

of the soil surface with the overlying atmosphere. Ham et al. (1993) reported short-

wave optical properties of five plastic mulches used in commercial crop production

(Table 2-2). These properties exhibited a wide range and thus helped explain earlier

work by Waggoner et al. (1960) and Eker (1967). Black and clear plastic mulches

provide extremes in the absorbance of short-wave radiation. The five mulches

described in Table 2-2 were used under bedded conditions in a field and

temperatures were recorded at the 10-cm depth below the mulched surface in the

center of the bed (Ham et al., 1993).

Table 2-2. Short-wave optical properties of five plastic mulches used in commercial
crop production. Data adapted from Ham et al. (1993).

Mulch Reflectance Transmittance Absorbance
Black 0.03 0.01 0.96
Sunfilm 0.12 0.37 0.51
Al-mulchi 0.39 0.01 0.60
White-on-Black' 0.48 0.01 0.51
Clear 0.11 0.84 0.05

1Al-Mulch and White-on-Black were Black mulch in which the upper surface was
painted with aluminum and white paint, respectively.

Diurnal measurements of soil temperature beneath the 5 mulches and a bare

soil plot as recorded by Ham et al. (1993) for a summer day are given in Fig. 2-1. The

highest temperatures were observed beneath the black mulch, which strongly

absorbed shortwave radiation. Cooler temperatures were observed beneath the

Al-mulch and White-on-Black mulches, which had lower absorbances but higher

shortwave reflectances. Although clear mulches are normally used to maximize soil

heating during cooler seasons, in this experiment the investigators placed the

plastic covers in intimate contact with the soil to ensure a minimal average air gap

between the mulch and the soil. The small air gap greatly minimized thermal

contact resistance at this location. Thus, the thermal insulation due to air between

the plastic sheet and the soil surface was minimized in that investigation.

38 Al-Mulch
U W White
34 Bare Soil


(; 26 -
Day of Year-195
22 0 4 8 12 16 20 24
Local Standard Time (hour)

Figure 2-1. Measurements of soil temperature at the 10-cm depth beneath five
different plastic mulches and a bare plot near Manhattan, KS on
day 195 of the year 1991. Adapted from Ham et al. (1993).

Besides mulching, the geometry of the soil bed, the directional orientation of

the beds, crop geometry (canopy structure) and the directional orientation of crop

rows also influence the radiation-energy exchange with the soil surface. The

amount of solar radiation reaching a horizontal unit of the surface of the earth (at a

particular latitude, longitude, and elevation) depends upon a number of factors

including the intensity of radiation emitted by the sun, astronomical considerations

determining the position of the sun, and the transparency of the atmosphere

F...i -r.-trir.: et al., 1983). The sum of direct solar radiation and diffuse sky radiation

received by a unit horizontal surface is referred to as global radiation. In humid

climates, clouds also play a major part in decreasing the amount of solar radiation

received at a location. The vegetative canopy of growing annual crops tends to

provide increasing shade for the soil surface as the growing season progresses, due


to interception of incoming solar radiation. The shadow of the growing crop results

in variable radiation interception at the soil surface, thus affecting the thermal

regime in the upper soil layers.

Coupled Heat and Water Flows

The energy state of water molecules in soil changes when a temperature

gradient is imposed on the system, and water moves in response to this energy

gradient. For the duration of the imposed temperature gradient, an irreversible and

spontaneous flow of heat occurs. Thus, heat and water flow processes interact with

each other and require mutual coupling. Where simultaneous movement of more

than one component (like water and heat) occurs in the system, isothermal flow

equations and simple equilibrium thermodynamics cannot be used to predict or

describe the water flow (Enfield, 1970).

If a soil system initially at equilibrium (uniform water content with a

constant water potential at constant temperature throughout) is subjected to a

temperature gradient, water will flow from a warm region to a cooler region (Gur et

al., 1952; Taylor and Cavazza, 1954; Woodside and Cliffe, 1959). When the degree of

water saturation in the soil is relatively high, water flow occurs primarily in the

liquid phase (Taylor and Cary, 1960). When the degree of water saturation is very

low, water movement occurs primarily in the vapor phase. At intermediate water

contents, however, movement occurs by a combination of liquid and vapor flows

(Cary and Taylor, 1962a, 1962b).

Possible magnitudes of thermally-driven flows have been summarized by

Cary (1966) for various soils and commercial materials. Vapor flow in soils in

response to thermal gradients of 0.01 C were observed to vary from about 0.4 to

2.0-mm day', the same flow rates that would require pressure-head gradients of

0.05 to 2.50 m of water. Since vapor flow out of wet soils (evaporation) may range

from near zero to perhaps 12 mm day' depending on climatic conditions, one can

see that temperature-induced vapor flow should not be ignored.

Philip-de Vries Theory (PDT)

Approximately 40 years ago, in an attempt to describe the phenomenon of

simultaneous flow where water is moved by a thermal gradient, Philip and de Vries

(1957) and de Vnes (1958) proposed theory and equations to describe coupled heat

and water transfer in porous media, in terms of classical mechanisms of vapor

diffusion and liquid movement by capillarity accounting for heat flow and its

resultant effects. Philip and de Vries (1957) were apparently the first to include the

interaction of vapor, liquid, and solid phases and the difference between average

temperature gradients in the air-filled pores and in the soil as a whole. Taking these

factors into account, an approximate analysis was developed which would predict

orders of magnitude and general behavior that was in satisfactory agreement with

the experimental results of Gurr et al. (1952), Jones and Kohnke (1952), Rollins et al.

(1954), Smith (1943), and Staple and Lehane (1954).

Philip and de Vries (1957) separated the isothermal and thermal components

of vapor transfer to show the effect of relative humidity. They separated the vapor


flux into two components: 1) one due to temperature gradients; and 2) another due

to moisture gradients. In the saturation range where liquid transfer occurs, the

liquid flux for water was separated into components caused by temperature

gradients, water gradients and gravity. Philip and de Vries (1957) suggested that

interaction between vapor and liquid phases occurs as a series of evaporation and

condensation steps coupled with liquid flow through so-called 'liquid islands'. When

a -ei .,ei.arur aiir ir .- i pi.-1 across a soil, water evaporates at the warm end of

the pores and diffuses as vapor, due to a vapor density gradient, towards the cooler

end where condensation occurs. This process changes the curvature of the liquid

menisci of a liquid island, which continues until capillary liquid flows through the

island, p:..,.:-d b. rhe ,-i.:, in.r difference of both curvatures, equals the rate of

condensation and evaporation. The overall effect is a decrease in the tortuous path

length and an increase in the cross-sectional area available for vapor diffusion.

Philip and de Vries (1957) presented the simultaneous equations of

continuity describing water and heat transfer under the combined influence of

water and temperature gradients as

S=V(DOVT) +v(DV6e) Kk aT = V(XVT) -Lv(Dve) [2.8]
t T -at

where Dr" is thermal water diffusivity, Do" is isothermal water diffusivity, k is the

vertical unit vector, Do' is isothermal vapor diffusivity, and L is the latent heat of

vaporization. Note that the existence of soil-water hysteresis may, however,

invalidate the use of Eq. [2.8] (de Vries, 1958).

Application of Philip-de Vries Theory

Various applications and verifications of PDT have been documented. Barnes

and Allison (1984) applied PDT for non-isothermal conditions to develop a model for

distribution of deuterium and 80 in dry soils. The model was later verified by the

same investigators (Barnes et al., 1989), with the experimental data showing

satisfactory results. However, they emphasized that the assumption of constant

tortuosity for both liquid and vapor movement is questionable. They reported that, m

moderately dry soils the liquid tortuosity factor is probably reduced, and vapor

diffusion may be enhanced by the presence of 'liquid islands'. Scanlon (1992) used a

model developed with PDT to evaluate liquid and vapor water flow in desert soils of

Texas using 3"C1 and 3H tracers. Observed and predicted potentials with the model

showed good agreement. Scanlon (1992) concluded from simulation results that

thermal vapor fluxes were downward-directed during the day and upward at night

during both winter and summer. Net vapor flux was also downward during the day

and upward at night.

Bach (1992) conducted isothermal and non-isothermal laboratory

experiments to investigate the significance of non-isothermal water flow and to

examine theoretical and numerical descriptions of the transport processes. He

concluded from sensitivity analysis that, with adjustment in Ch (temperature

coefficient of the matric potential), PDT provided an adequate description of the

non-isothermal transport processes.

Many other applications of PDT to non-isothermal water movement have

been observed to describe and simulate transport processes (Benjamin, 1989; Chung

and Horton, 1987; de Silans et al., 1989; Mahrer et al., 1984; Milly, 1984; Nassar and

Horton, 1992; Noborio, 1995; Sophocleous, 1979; and Yakirevich et al., 1997). Cassel

et al. (1969) concluded that PDT showed acceptable agreement with observed soil

water movement.

In order to account for soil inhomogeneity and hysteresis of the moisture

retention process, Milly and Eagleson (1980) modified the PDT by using soil matric

pressure-head gradients instead of soil moisture gradients. These authors also

included effects of the heat of wetting on transport processes. A finite-element

one-dimensional formulation of the governing equations showed that the model

performed well in simulating highly coupled, hysteresis-affected, or very nonlinear

problems. Later, Milly (1984) used this modified form to study soil-water evaporation

from silt loam and sandy soils and reported that soil water evaporation was

generally more sensitive to isothermal than to thermal vapor diffusion.

Published evidence thus shows that the mechanistic PDT appears to be

accurate for coupled heat and water transfer in soils. This mechanistic approach to

formulate expressions for the fluxes of water and heat is based upon concepts

derived from fluid mechanics, mass diffusion and heat conduction.

Solute Transport during Coupled Heat and Water Flow

Simultaneous flow of solutes and water occur in most soil-water systems

under field conditions. Salt concentration gradients may also be a major contributor

to water movement in soils located in ard climates. Salt concentration gradients

often occur within the uppermost soil profile, where evaporation selectively removes

water at the soil surface and result in an accumulation of salt. Kemper and Rollins

(1966) reported that, for experiments involving Wyoming Bentonite clay, soil

moisture suction heads needed to be greater than 50 m of water before the osmotic

potential gradient would c ,.- ,1iarfi,.: in r:. ii water movement. Low (1955)

observed that, when the solute was completely restricted by the porous medium,

osmotic potential differences were as effective as pressure potential differences in

causing water to move through clay under isothermal conditions. In studying water

flow through a porous ceramic filter, Jackson (1967) reported that soil water flow

could deviate from Darcy-type flow due to the presence of salt. In laboratory studies

on clay loam and sandy loam soils, Qayyum and Kemper (1962) reported that

diffusive movement of water occurred to the salt-bearing surface at soil moisture

contents less than one-half of field capacity, and that viscous flow occurred at soil

moisture contents higher than one-half of field capacity in closed soil columns.

Kemper and Maasland (1964) found that salt sieving increased with

decreasing pore size of the soil and decreasing valence of the saturating cation.

Letey and Kemper (1969) introduced a theory to describe simultaneous movement of


water and solute in soil. The theory included an equation to describe water flow and

another equation to describe solute transfer. The equation for solute transfer implied

that the movement of solute is controlled by Fickian diffusion and salt sieving.

In summary, temperature and osmotic gradients can induce water flow

within soils. As water moves, heat and solute are transported, thus altering the

driving forces (gradients). Interactions among thermal, soil moisture, and osmotic

effects exist, and they differ for the liquid and vapor phases (Nassar and Horton,


Nassar and Horton (1989b) extended the PDT, describing water transport

under the combined effects of matric and osmotic pressure head gradients in

unsaturated, non-isothermal, saline soil. Their theory was based upon the concept

of viscous flow of liquid water under the influence of capillary, adsorptive, and

osmotic forces and on a concept of vapor movement by diffusion. They presented

two equations for vapor and liquid transfer as

S-D= VT DevL + DcvVc [2.9]
S=-D VT DVeL + DcVC Kk [2.10]

where q, is water vapor flux, D, is thermal vapor diffusivity, DTL is thermal liquid

li rj I -, r.. Dev is isothermal vapor diffusivity, Def is isothermal liquid diffusivity, Dcv

is vapor diffusivity due to solute concentration, DcL is liquid diffusivity due to solute

concentration, and c is solute concentration.


Equations [2.9] and [2.10] were summed in order to describe the total water

flux as

q = -DVT DeV + DVC Kk [2.11]
where q, = q, + q, is total water flux, DT" = D, + D, is thermal moisture

i ffu,; ir r, -,w" = Dev + DOe is isothermal moisture diffusivity, D. = Dcy + DcL is

moisture diffusivity due to solute concentration, and 6 = O, + O6 is total volumetric

water content, where vis water vapor content (expressed as equivalent liquid


Applying the principle of heat conservation, an equation for transient heat

transfer combining the effects of matric pressure head, osmotic pressure head, and

temperature was formulated. For developing the solute-transport equation (non-

interacting solute transfer), Nassar and Horton (1989b) included the mechanisms of

molecular diffusion, hydrodynamic dispersion, thermal diffusion, mass flow of water

(convection), and salt sieving.

Nassar and Horton (1989a, 1989b) conducted experimental and numerical

studies of coupled heat and mass transport in closed soil columns with three

different soil textures. Experimental observations for unsaturated conditions showed

water accumulation within the low-temperature region of the columns, but solute

accumulation within the high-temperature region. Predicted and observed results

were in good agreement when initial soil-water contents were relatively large and

initial solute concentrations were low. However, for small initial soil-water contents,

the model underestimated soil-water contents close to the low-temperature

boundaries and overestimated than near the high-temperature boundaries. For high


initial solute concentration, the theory underestimated liquid transfer toward the hot

boundary temperatures. Unfortunately, the theory proposed by Nassar and Horton

(1989b, 1992) is restricted to homogenous soils and cannot be used in its presented

form for heterogenous soils. The water flow equation is written in terms of water

content and thus does not allow simulation for multilayered soil. Furthermore, the

model is primarily restricted to soil conditions of water unsaturation.

Temperature Effects on Solute Transport

The rates of most chemical reactions increase with temperature, changing

their equilibrium constants (Clark, 1996). Clark further states as a rule of thumb,

that the reaction rate doubles for each increase in temperature of 10 C. Generally,

the understanding of a chemical reaction is not considered complete until the

sensitivity of the reaction to temperature is understood. Many chemical studies thus

include experiments at different temperatures.

Daniels and Alberty (1961) and Clark (1996) discussed the influence of

temperature on reaction rate in the following way. Kinetic data over a range of

temperatures can usually be represented by use of an empirical equation first

proposed by Arrhenius
kr= F e RT [2.12]

where k, is reaction rate, F, is a frequency factor, E, is activation energy, and R, is

the universal gas constant. Set in logarithmic form,

logk, + log Ff [2.13]

According to Eq. [2.13], a straight line should be obtained when logarithm of the

rate coefficient (k,) is plotted against the reciprocal of the absolute temperature (T).

The E, can then be estimated as -2.303 R, x the slope of the line. The activation

energy may be envisioned as an energy barrier that must be surmounted for the

reaction to proceed. An increase in E, indicates a slowing down of the reaction (and

vice-versa). The E, thus can be determined by measuring the reaction rate at

different temperatures.

The influence of temperature on chemical processes has also been described

by Stumm and Morgan (1981) through use of the Arrhenius equation. Simfnek and

Suarez (1993) described a function with some modification to account for relative

change in rate coefficients with temperature as

m r re J [2.14]
f(T) = e '[2.14]

where T, is the reference absolute temperature and f (T) = 1 at T = T,.

Cohen and Ryan (1989) and Cohen et al. (1988) used the relationship given by

Bird et al. (1960) to describe the temperature dependence of diffusion coefficients.

The chemical molecular diffusion in air as a function of temperature was estimated


Da (T) = DT) 1.823 [2.15]

where D,(T,) is the chemical molecular diffusion coefficient in air at the reference

temperature. Temperature variation of the chemical molecular diffusivity in liquid

water was estimated as

D ( T) = D T T) [2.16]

where p, is the viscosity of liquid water.


The PDT was proposed !:1, aj- ij:o to mathematically describe coupled

water and heat flows in soil. It has been well-verified and applied by a number of

researchers to date. The theory was extended by Milly and Eagleson (1980) to make

it applicable to heterogenous media, by replacing the water content-based

formulation with a matric potential-based formulation of flow equations. This

advancement was a big step, as it made the theory suitable for the solution of field

problems. Nassar and Horton (1989a, 1992) extended the theory to include salt

effects on water flow and provided insight into additional coupling. But, this theory

cannot be used to solve field problems, because Nassar and Horton used the

original PDT (water content-based) formulation for modification.

Relatively little evidence exists that the potential of PDT and its modification

by Milly and Eagleson (1980) has been fully exploited to use the improved estimates

of water content/flow, and temperature information. Several researchers (Barnes

and Allison, 1984; Barnes et al., 1989; Nassar and Horton, 1989a, 1989b; Noborio,

1995; Scanlon, 1992; and Yakirevich et al., 1997), however, utilized PDT by making

use of improved estimates of water content/flow when exploring solute transport.


These investigations, however, did not attempt to use the improved estimates of

temperature by inclusion of its effect on solute transport.

A need, therefore, exists to exploit the full potential of PDT in solving field

problems where temperature could be highly transitional both in space and time.

Furthermore, where chemical (especially volatile-compound) transport is also of

concern, a need exists to use both the improved estimates of water content/flow and

of temperature provided by PDT. The work reported here utilizes PDT in totality to

develop a 2-D numerical model for solving problems involving the fate and transport

of volatile chemicals that may exist in all three phases (viz. liquid, solid, and gas) of

soils. Details for the developed model are presented in Chapter-3.


The model developed here considers the soil to be a heterogenous isotropic

medium composed of solid, liquid, and gaseous phases. The latter is assumed to be

a non-interacting binary mixture of water vapor, dry air, and volatile chemicals.

Other general assumptions in the model include

1. The solid matrix is undeformable and immobile.

2. The: i-,jii ... ._er i r ILr: ':. i lP, LIri

3. The density and viscosity of liquid water are not influenced by dissolved


4. Hysteresis in the relationship between capillary-pressure and hydraulic

conductivity is neglected.

5. Liquid water and water vapor are in local thermodynamic equilibrium.

6. Water vapor and other gases behave as perfect gases and advective transport

of these components is neglected.

7. Local thermal equilibrium exists among all phases.

8. Radiative heat transfer is negligible.

9. Thermal and matric potential gradients interact to provide heat and water

flow in the porous medium.

More specific assumptions are stated as required. Additional details

concerning the equations governing flow and transport are presented in the


Liquid and Vapor Water Flow in Soil

The original water content-based formulation of PDT was modified by Milly

and Eagleson (1980) to provide

[1 P + L eaaPj ah [ +D,,)Vh +D,VT+Kk

pj + e, aPI aT", a Tv
1- +--
p, 9T p, 9T 9t

where p, is vapor water density, 0O is the volumetric air content, and D, is

isothermal vapor diffusivity. Equation [3.1] describes both liquid water and vapor

water flow in response to gradients of matric potential and thermal energy.

Heat Flow in Soil

Coupled heat flow in the soil in response to gradients of temperature and

pressure-head has been described by (de Vries, 1958; Milly and Eagleson, 1980)

C+H +p H aT = VT+p LD h C(T-T q [3.2
S(T 2 T ) t L hO J [3.2


C =Cd+CLpL (I+Cvp HI = [Lo+C,(T- To

H2 (CLpL \Cp(T T) p Lo [3.3]

where Cd is volumetric dry heat capacity of the total medium; CL and C, are specific

heat of liquid water and vapor water; respectively; and the subscript implies the

magnitude at a designated reference level.

Transport of Multiple Species of Solutes in Soil

Following the approach of Similnek et al. (1992) and Simfnek and van

Genuchten (1994), in which solutes can exist in all three phases (liquid, solid, and

gaseous) with non-equilibrium transport, -i r.:.- -ir. rg transport equations can be

stated as

a + +-c= + aDELc c 6aDlgg c
a a t -at vLDvc+v Dl g-V(qLc1l) [

+E)LC= + =.) +) D VC DVg -VL i -pL, ri + C -

s, + Ps, ,i p+s 1 ie g ,i-l i1 s, i-lpsi-1 [3.5]

Plg, i.i-1i+YL, iL+s, iP +g, a; i=2,..., n (species)

where s represents adsorbed concentrations; g is the gas concentration; PL, P., and

/, are the first-order rate constants for solutes in the liquid, solid, and gas phases,

respectively; YL, y,, and 7, are zero-order rate constants for the liquid, solid, and gas


phases, respectively; p is the soil bulk density; D1 is the dispersion coefficient for the

liquid phase; and D9 is the diffusion coefficient for the gas phase. The subscnpts L,,,

and correspond with the liquid, solid, and gas phases, respectively; the subscript,

represents the i" chain number; n, is the number of solutes involved in a chain

reaction; and p,*, p,*, and p,* are first-order rate constants for solutes in the liquid,

solid, and gas phases, respectively, providing the connection between individual

chain species. The nine zero- and first-order rate constants in Eqs. [3.4] and [3.5] can

be used to represent a variety of reactions or transformations including

biodegradation, volatilization, precipitation, etc., and the effect of temperature on

these coefficients is included through use of Eq. [2.14]. E-u r i.:ri [3.4] alone is

sufficient for a single solute that does not undergo transformation to related

chemical species in a 'chain' reaction. ?..rh E 1- [3.4] and [3.5] are required when the

primary solute undergoes transformations to other species or when multiple

unrelated solutes are present in the soil system.

In order to solve the solute-transport equation, a knowledge of liquid water

content, air content, and volumetric flux is required. Solution of the coupled water

and heat flow equations (Eqs. [3.1] and [3.2]) provides this information.

Enerqv Balance at the Soil/Atmosphere Interface

Interception of solar radiation, interaction with the overlying atmosphere,

and loss of water through evaporation is described via an energy-balance equation

for the soil surface


(1-PR +R RLA = H+LE+G [3.61

where p, is soil surface lib-.J i, measured global radiation (taken as positive

downward), R, is longwave sky irradiance (positive downward), Rs is longwave

radiation from the soil surface (positive upward), H, is sensible air heat flux (positive

upward), LE is latent heat of the evaporative flux (positive upward), and G is soil

heat flux (positive downward).

The solution of Eq. [3.6] provides information for soil surface temperature and

atmospheric evaporative demand (excluding transpiration) required to specify the

top boundary condition of the solution domain. This equation is strictly valid only for

a soil surface without a mulch cover.

The presence of a mulch cover alters interaction of the soil surface with the

overlying atmosphere (as discussed earlier). An approach based on the work of

Hares and Novak (1992) and Ham and Kluitenberg (1994) was adapted for

conditions where the soil surface is covered with a plastic mulch such as commonly

used in soil-bed systems for vegetable production in Florida. The energy balance for

a plastic mulch can be defined as

Rnm + Hs Cc = 0 [3.7]

where Rn, is net radiation of the mulch surface and C, represents conduction or free

convection between the mulch and free soil surface. Energy balance for the soil

surface beneath the thin plastic mulch can be defined as

Rn + C+ G = 0 [3.8]

where Rn, is net radiation of the soil surface. These energy-balance equations

neglect heat storage within the thin plastic mulch, and the potential impacts of cuts

and planting holes m the mulch is not considered. Thickness for such plastic

mulches typically ranges from 0.02544 to 0.15264 mm.

Net radiation of the plastic mulch is given by (Ham and Kluitenberg, 1994)

Rn. = aRg 1+ P' P P r + 1 ) + Pr msTs

1-e) -2eT, [3.9]
where T,, T,, and T, are the temperatures of the sky (usually taken to be

atmospheric temperature at 1 m above the ground surface), soil surface, and mulch,

respectively; e,,y, e,, and em are the emissivities of the sky, soil surface, and mulch,

respectively; a. is the shortwave absorbance of the mulch; T,, is the transmittance

of the mulch in the longwave spectrum; and a is the Stefan-Boltzmann constant.

The variables p* and p,* represent the internal reflection functions for shortwave

and longwave radiation, respectively. These functions account for multiple

reflections of radiation between the mulch and soil surface, and are defined by Van

de Hulst (1980) as

P' = 1 P ) = 1 + P P p [ 3.10]

S 1- 1- = p, 1- + 2 [3.11]

where p, is the shortwave reflectance of the lower mulch surface and p,, is the

reflectance of the mulch in the infrared spectrum. The infinite series expansion

shows here that the reflection functions account for the initial radiation in the layer

and then amplify the effects via repeated reflections between the soil surface and

the underside of the mulch.

Sensible heat transfer between the plastic mulch and the atmosphere is a

function of the aerodynamic properties of the surface and the temperature gradient

between the plastic and the boundary layer, as defined by

C ( T T )
Hs = C(T- [3.12]
a, h

where Ta is air temperature and r,,h is the aerodynamic resistance to heat transfer.

Heat transfer between the mulch and the underlying soil surface is estimated using

the equation
C [3.13]

where r, represents a thermal contact resistance at the mulch-soil interface and

mainly depends on the air-gap thickness between the soil surface and the mulch.

This determines whether the heat transfer is conductive or convective. Ham and

Kluitenberg (1994) described the estimation procedure and presented representative

values of r. for different air gaps between the soil surface and the mulch for various

plastic mulches.

Net radiation of the soil is dependent upon the global and longwave sky

irradiance transmitted through the mulch, radiation emitted from the mulch and

soil, and multiple scattering of radiation between the soil and mulch as given by

(Ham and Kluitenberg, 1994)

n,= (-I P)TP'Rs pre (Im4 a,,ersoT kPY,,i T ,)-eT [314]

Given the subsurface temperature near the surface at depth z, ground heat flux can

be approximated at the future time step p + 1 (p being the current time) as (Horton

and Chung, 1991)

G=- T' f Ts-- T1 + C s Ts ]z [3.15]
Az s 2At

where T, is soil temperature at the first node below the soil surface and (Tf'i

T,)IAt is the rate of change in temperature with time at the surface between time

steps p and p + 1.

Equations [3.1]-[3.5] require supplementation with appropriate initial (IC) and

boundary (BC) conditions to describe a specific simulation scenario. The IC in

general (at t = 0) are specified as h (x, z, 0) = h, (x, z), T (x, z, 0) = T,(x, z), and

c (x, z, 0) = c, (x, z) respectively for Eqs. [3.1], [3.2], and [3.4], [3.5]. Here h,, T,, and c,

are known values of the state variables at the beginning of the simulation; and x

and z represent the horizontal and vertical locations, respectively.

Boundary conditions may be defined as a specified value of matric

potential/or temperature and/or solute concentration (Dirichlet-type BC), a specified

water/or temperature and/or solute flux (Neumann-type BC), or a combination of the

two (Cauchy-type BC). The top boundary (i.e., the soil surface) in the system under

consideration is a special case of the Neumann-type BC for water (evaporative flux),

dependent upon the overlying atmosphere. In general, the boundary conditions

depend on the specific system being investigated.

Input Parameter Requirements

Soil input data include sand, silt, clay, and organic matter contents; bulk

density; saturated hydraulic conductivity; and soil hydraulic parameters which

describe the water retention curve (Eqs. [3.7]-[3.8]). For the soil/atmospheric

boundary conditions, weather data such as average air temperature, air temperature

daily amplitude, average dew-point temperature, dew-point temperature daily

amplitude, wind speed, surface roughness length, and day-length are required

inputs as well.

Soil hydraulic properties are estimated following van Genuchten (1980) as

e5- r 1
e = e + [ -(: ; FM = 1 -
0+ [ +(.:. I 1

K(L) Kat [3.16]


h 1 -) e -e [317 1

where 6, is residual water content, 0, is saturated volumetric water content, |h is

the absolute value of the matric potential, and K_, is the saturated hydraulic

conductivity. The parameters a, rL, and rr) are empirical constants used to fit the

moisture retention curve. The temperature dependence of matric potential (h) was

described as (Philip and de Vries, 1957)

dh h do,
dT = hdT Chh [3.18]
and hydraulic conductivity (Kh)) a(Constantz 1982)
and hydraulic conductivity (K(h)) as (Constantz, 1982)

K(h) Kef(h) [3.19]
w PL, ref

where o, is surface tension of water; Ch ( = -0.0017 K-' by Milly and Eagleson, 1980)

is temperature coefficient of matric potential; p,"' and p T are dynamic viscosities

at the reference temperature and soil temperature, respectively; pL,r and PL, are the

densities of water at reference temperature and soil temperature, respectively; and

K,, (h) is hydraulic conductivity at the reference temperature.

The dependence of 6L upon temperature can be determined as (Milly and

Eagleson, 1980)
ae, 8eL
S h exp[-C,.(T- To [3.20]

Water vapor density is given by Marshall and Holmes (1979) as

Pv = gaPLI u T [3.21]

where p is water vapor pressure and M, is the molecular weight of water. Water

vapor pressure is related to soil water matric potential by the relation (Marshall and

Holmes, 1979)

S= 0, exp [hM] [3.22]

where cp, is the saturated vapor pressure, given by Noggle (1985) as

S= 237.7 exp '~ -- [3.23]

Differentiating Eq. [3.21] with respect to h yields

ap M mJ [T .324
ah 9 ga- PaL uT exp gahRMT [3.24]
R wT

and with respect to T gives

ap, M
aT D1 *D2- g PL R T [3.25]

Dl D2= D3+qsah M D3 ; D3=exp^ gAh

Variation in saturated vapor pressure with temperature is given by

a0, 237.7ML I ML
aT 2 ex+p [3.26]

The vapor diffusion coefficients are determined as (Milly and Eagleson, 1980)

Dh ,DTe D T -- [3.27]

2 f = for 8L<60 f = 6(+ T- ( for ,L>e, [3.28]

where D, is a molecular diffusion coefficient of water vapor in air, T' is an air-phase

tortuosity factor, fis vapor path correction for liquid islands, 0T is total porosity, 8, is

the water content at which liquid continuity breaks during desaturation, and

S( = 1.7) is a temperature-gradient enhancement factor for vapor.

Specific heat is estimated from the methods of deVries (1963) as

C = X c +X c +Xcc+Xo com +CL + Cve [3.29]

where X,, X,, X, and X.m are the volume fractions of sand, silt, clay, and organic

matter, respectively; and c,, c, c0.,, C, are the volumetric heat capacities of sand,

silt, clay, organic matter, and water vapor, respectively. The specific heat of water

vapor is given by Noggle (1985) as

C 0.0180015 [3.30]

The apparent thermal conductivity can be estimated as a harmonic mean of

the soil constituents (de Vries, 1963)

El Xi. -
1.25 i=st,c,o, k 1+ -I-1
E k*Xi 3 -j= .,4 X gj
1=s,st,c,o,8 J
where X, represents the respective thermal conductivities of the soil constituents, g,

are shape factors for the soil constituents, and X, is the thermal conductivity of air.

The incoming global radiation for the soil/atmospheric boundary condition

must be partitioned into a flux-density distribution for daylight hours. Chung and

Horton (1987) suggested a sine wave distribution

R R sin t-SN+ DL3.32]
g 2 DL 2 DL

where Rd. is daily global radiation, DL is day length, SN is solar noon, and t is the

time of day. The air and dew-point temperature distributions can also be expressed

as sine functions

S ( 2nt __ 2nt
T = T,+A sin 86400 +n = Td+d Asin 4 +0 j [3.331

where the bar indicates average temperature; T. and Td are air and dew point

temperatures, respectively; A, and Ad are air and dew point temperature amphtudes,

respectively; tm is the time from midnight; and the maximum temperature is

assumed to occur at solar noon.


The soil surface emissivity, e,, and albedo, p,, are described by the method of

van Bavel and Hillel (1975, 1976) as

es = 0.9+0.186L p, = 0.25 for eL-0.10,
P, = 0.35-eL for O.lO-.e<0.25, pS = 0.10 for eLo0.25


Numerical solutions to differential equations require testing to ensure that

numerical instabilities have not been introduced into the solution due to round-off

or truncation errors and that the solution converges on a reasonable estimate of the

exact solution. A numerical scheme is stable if, with increasing time, the difference

between the numerical solution and the exact solution tends uniformly to zero at all

nodes (Benjamin, 1989). Convergence is satisfied if the numerical solution

approaches the exact solution as spatial- and temporal-step sizes approach zero.

Proper testing also ensures the correctness .f ri -ri, i-, ] c3L] r-~pe-e-r1r-1 ,: tr,-.: the

physical system being modeled and the validity of any assumptions made in the


Testing of a model may be achieved in multiple ways. One such way is by

comparing the model results with analytical (exact) solutions for greatly simplified

conditions, and another is by comparing the model results with observed data from

carefully designed experiments. Unfortunately, no analytical solution to the

governing equations for coupled heat and water transport is known to exist. Hence,

portions of the model were tested with observed data and analytical solutions for

simplified conditions, where available.

Isothermal Water Flow

Model results were compared with data for transient one-dimensional,

isothermal water flow through sand as reported by Haverkamp et al. (1977), with

vapor flow ignored. Haverkamp et al. (1977) reported K., 6,, 6,, and soil bulk

density to be 34.0 (cm h '), 0.287, 0.075, and 1.66 g cm-3, respectively. An analytical

model by van Genuchten (1980) was fitted to the soil water content versus matric

potential data in order to provide optimum values for the hydraulic parameters a

and r. A uniform initial water content of 0.10 was specified and water content at

the 0.70 m depth was kept constant at this value. A constant water influx of

13.69 cm h' was specified at the top boundary. Results of the numerical simulation,

and experimental data for different times, are shown in Fig. 4-1. The numerical

solution was stable and converged upon a close approximation of the experimental

values. The model accurately described the movement of the water front in the

column of sand.

,0 5 0.10 0.15 0.20 0.25 0.3
o observed.
-10 Estimated -1o
o 0
-20 0.2 hr -20

-30 -30 Figure 4-1. Observed and computed
S 0.4 r distributions of water content in a
-0 : -40 sand column during constant flux
S-50 -50 infiltration.

Water Content (cm3 cm3)

Uncoupled Heat Flow

Next, the numerical solution for heat flow was verified with field-measured

temperature data (Horton and Chung, 1991) for a Nicollet clay loam soil profile with

a bare surface located near Ames, IA. The effect of water flow was ignored.

Representative values of heat capacity and thermal conductivity for loam soil were

used. Initial conditions as specified by Horton and Chung (1991) were used in the

model and temperature at the 60-cm depth was held constant at 20 "C.

Temperatures observed for the soil surface were used in the model to provide the

surface boundary condition. One-dimensional numerical results (Fig. 4-2) for the

5- and 15-cm soil depths provided satisfactory representation of temporal

temperature profiles for each depth.

40, 20 40 60 80 100 12o
Observed 5-cm
36 Estimated 5-cm 36
o Observed 15-cm
U 32 ------ Estimated 15-cm 32

28 28

24 24
E : /<
20 o -20

16 16
0 20 40 60 80 100 120
Time (hours)

Figure 4-2. Observed and computed temperatures for two depths
in a field soil.

Coupled Water and Heat Flow

In order to test the capacity of the numerical solution to predict coupled heat

and water transfer, simulation using the experimental data of Nassar and Horton

(1989a) was performed. Nassar and Horton (1989a) conducted laboratory

experiments of heat and mass transfer with controlled temperature gradients m

sealed columns (14 cm in length) of silt loam soil with initially uniform water

content (0.117). The column ends were subjected to different temperatures (19.03 C

and 9.36 C) in order to create a temperature gradient. The experiment was run for a

period of 31 days. The observed mean value of saturated hydraulic conductivity for

the soil was 0.954 cm hr-', with a corresponding soil water retention curve having

been determined by Nassar and Horton (1989b).

The van Genuchten (1980) model was fitted to the soil water retention curve

in order to obtain optimum hydraulic parameters for the soil (a = 0.03137 cm',

Ir = 1.33308, 68 = 0.65, 6, = 0.05) for silt loam soil. A numerical simulation was

carried out for 31 days using zero flux boundaries at both ends for water flow and

Dirichlet boundaries for heat flow. Observed and predicted water contents are given

in Fig. 4-3. Good agreement was noted everywhere except in the vicinity of the cold

end of the soil column with elevated water content. The predicted total water

content in the column was slightly higher than the observed value. Since the

computed mass balance error was less than 1.0 % in the numerical simulation, the

possibility exists for a small systematic error in the measurement of either initial or

final water contents during the actual experiment.

------- Estimated
0.25 Observed
S-Initial Figure 4-3. Observed
E0.20 and computed water
Content after 31 days
S0.15 under the influence
\ of both matric head
00 and thermal
) 0.10
S01 0o gradients.

000 2 4 6 8 10 12 14
Distance (cm)

Isotherm al 1-D 1 -,Iur- Tr .r..... ri. unri _rea.1, ''J r .r l.:....

The model component for solute transport was tested against data reported

by van Genuchten (1981) for the movement of a pulse of boron 'H ,E,-', during steady

water flow through a column of Glendale clay loam. The solute-transport

parameters (imfmnek and van Genuchten, 1994) used were: DL = 2.042 cm2 hr1,

68 = 0.445, p = 1.222 g cm-3, kd = 1.14 cm3 g', P = 1.0, T1 = 0.0, f, = 0.47, and

o) = 1.33x10-2. The initial concentration was zero and a pulse of aqueous solution

containing a boron concentration of 20.0 mmol, L' was applied at a constant flux of

713.0 cm hr-' for a duration of 121.44 hr. Assuming two-site chemical non-

equilibrium sorption, a simulation was run for 20 pore volumes of effluent. A

comparison of numerical and experimental results (Fig. 4-4) verified that the model

described the transport process with satisfactory accuracy.





Figure 4-4. Observed and computed break through curves for boron.

Isothermal 2-D Solute Transport during Steady Water Flow

The mathematical accuracy of 2-D numerical simulations for isothermal

steady water flow, was verified employing a simple analytical solution. Solute (non-

volatile) transport in a homogenous, isotropic porous medium during steady-state

unidirectional groundwater flow can be written as

ac a2c 2c ac
R = D D R c 4.1]
at TX2 LZ2 az Rc

where D, and DL are transverse and longitudinal dispersion coefficients,

respectively, and v, is average pore-water velocity in the downward vertical

direction. The initially solute-free medium is subjected to a solute source, C,, of unit

concentration. The source covers a length 2a along the inlet boundary at z = 0, and

is located symmetrically about the coordinate x = 0. The transport region of interest

is the half-plane (z 0 0; -- < x s o-). The boundary conditions can be written as

c(x,0,t) = Co -asx a; c(x,0O,t) = 0 other values of x

i, mc ac [4.2]
lim- = 0 ; 1. 0
z-_ az .. .. .

An analytical solution to the above transport problem is given by Leij and Bradford

(1994) as

Cz ( O 2 z2
c(x,z,t) Cz -- exp iR+ z
4(xDL), 2DLJ [ 4D, 4DL0,rW
a -x er a+x jid
2(DT } )l [ 2(Dr.w)a)"

0 _3 ___ ___ _
SOURCE (a =50-m)
-10 Analytical
-20 -50 Days

-30 1 ___ _0 ays
2 Days
-40 lConcentration contour level= 0.1
Sv= 0.1 m day'

D = 1.0 m2 day'
R= 1.0
-60- C= 1.0
0 20 40 60 80
Distance from mid-source (m)
Figure 4-5. Location of solute concentration fronts for three times as estimated
numerically and with the analytical solution.

The numerical simulation assumed a source-width of 2a = 100 m. Because of

symmetry, computations were carried out only for the half plane where x 2 0.

Figure 4-5 shows the estimated solute concentration fronts (0.1) at selected times


and various values of the parameters as used in the computations. F i- Il p rre.1'- r-.

by the numerical model were similar to those given by the analytical solution.

Isothermal 1-D Solute Transport with Nitrification Chain during Steady Water Flow

To verify the numerical solution capacity in simulating multiple linked

species in a chain, comparison was done with results generated from an analytical

solution published by van Genuchten (1985) for 1-D transport of solutes involved in

sequential first-order decay reactions. The analytical solution holds for a

homogenous, isotropic porous medium with steady-state uni-directional

groundwater flow. The solute-transport Eqs. [3.4] and [3.5] for this case reduces to
ac, = a2c ac
R- D vx PL,1 R cl [4.4]
R1 at ax2I44]
ac. 2 82c ac
R = a2 Vx aX PL,-1Ri-i- -1c PL,iRici ; i = 2,3 [4.5]

A simulation was performed here that involves a three-species nitrification chain

NH' (i=1) NO2(i=2) NO, (i=3) [4.6]

and follows van Genuchten (1985) and Cho (1971).

A NH,4 solution was applied to an initially solute-free medium. The input

transport parameters used for the simulation were

v = 1.00 m h' R, = 2.0 [-]

DL= 0.18 m2 h' R = 1.0 [-]

S= 0.01 h- R, = 1.0 [-1

2 = 0.10 h- c, = 0.0 [-]

P3 = 0.0 h-' o, = 1.0 [-]

The simulation was conducted for 200 hours. Figure 4-6 shows the

concentration profiles for all three solutes at time 200 hr, calculated both

numerically and with the analytical solution from the CHAIN code of van

Genuchten (1985). The numerical solution duplicated the analytical results almost


1.0 t------------------
-- Numerical NH*
o Analytical NW Figure 4-6.
S0.8 -- Numerical NO;
Analytical NO; Analytically
Numerical NO; determined and
O Analytical NO; numerically
U computed
0.4 concentration
profiles at 200
0.2 hours.

0 40 80 120 160 200
Distance (m)

Isothermal 1-D Solute Transport involving Gas Flow

The capacity of the numerical solution in simulating gas transport was

tested against an analytical solution given by Crank (1975). Considering uniform

constant water content in a homogenous 1-D soil column, for a chemical species

existing in the gaseous phase undergoing hydrolysis (represented as first-order loss,

pg), the transport equation (Eq. [3.4]) reduces to

ag D -pg 47
Ot ax2


For an application of finite source (g,) at one end with zero flux BC and zero

gradient BC at the other end, an analytical solution is (Crank, 1975)

gxt) = e go erf- +erf 2/--- [4.8]
2 n=- 2Dt 2Dt

where a is the initial source length. Prior to deriving the inril,'t~.: ]i ii ir:in the

transformation g = g* exp(-pg t) was used to simplify Eq. [4.7] to

=t D 2 [4.9]

A numerical simulation was conducted for 5 hours with g, = 50 moles,

a = 1 cm, D9 = 47.52 cm2 hr-1, and P = 0.055 hr The numerical solution

duplicated the results of analytical solution at 1, 3, and 5 hours (Fig. 4-7).

Figure 4-7. Analytically determined and numerically computed concentration
profiles at 1, 3, and 5 hours.


All major components of the numerical model were tested for mathematical

accuracy, the adapted numerical scheme, and its ability to represent the physical

systems. Both analytical solutions for simplified cases and published data were

used to check against simulation results from the numerical model. Excellent

agreement was observed in nearly all the cases undertaken for testing of heat-flow,

water-flow, and solute-transport model components. The coupled water-heat flow

component was also successfully tested. The model could not be tested in its

entirety as a single unit, because no analytical solution exists at present for coupled

water-heat and solute transport. Also, no published data set was found to test all

the features of the model together as a single unit. However, successful testing of

individual model components provides confidence in the application of the model

simulations conducted in this work.


a n1-~rij Ti j-,i: L :.r oif T.I,-ih i r.n.,ij- FunuQ ir.,,r ,rl :,,, i, i

Methyl bromide (MBr) is an odorless, colorless toxic gas used in agriculture

throughout the world. The United States uses it on more than 100 crops as a soil

fumigant, a post-harvest treatment, and as a plant quarantine treatment to control a

variety of pests. MBr is also used to fumigate structures such as grain storage

warehouses, flour mills, and ships carrying agricultural commodities. California and

Florida are the primary users of MBr, followed by Hawaii.

MBr, a broad-spectrum pesticide, is urli z-.1 rI-' "rd, as a very effective

agricultural fumigant during vegetable and fruit production in Florida, California,

and Hawaii. MBr is very effective for eradicating insects, nematodes, fungi, bacteria,

and weeds. In the United States, strawberries (16 % of U.S. total) and tomatoes

(24 % of U.S. total) are the crops which use the most MBr, consuming about 6,500

tons annually. Other crops which use this pesticide as a soil fumigant include

tobacco, peppers, grapes, and nut crops (USEPA, 1997).

In the United States, the U.S. Clean Air Act requires that the production and

importation of any substance defined as an ozone depleter by the Montreal Protocol

of 1991 be phased out (USDA ARS, 1997). Therefore, despite its economic value, the


U.S. Environmental Protection Agency (USEPA) plans to phase out the use of MBr

by January 1, 2001. Since the Clean Air Act allows no exceptions, all uses of MBr,

including plant quarantine, will be banned. In addition to being a toxic gas, MBr is a

-..3r i i.:a'nt .:-z:, r,-depleting substance, with recent scientific evidence estimating

that bromine from this material is at least 50 times more effective at destroying

ozone than chlorine from carboflurocarbons on a per-molecule basis (USEPA, 1997).

MBr has an average residence time of 1.5-2.0 years in the atmosphere (Hacherl,

1994). However, there are also other natural sources, primarily forest fires and

marine organisms besides anthropogenic sources, which contribute significant

quantities of MBr to the atmosphere.

Table 5-1. Earlier estimates for source-contributions of MBr to the atmosphere
(Grojesan, 1991; Khalil et al., 1993; Lobert et al., 1990; Mano and
Andreae, 1994; Yagi et al., 1993; and Zurer, 1993).

MBr Source Quantity in atmosphere
(m ton y )
Oceanic 30,000 (30 %)
Bio-mass burmng 30,000 (30 %)
Industrial 15,000 (15 %)
Agricultural 25,000 (25 %)
Total 100,000

Quantification of MBr release to the atmosphere is extremely difficult, since

limited data exist for the total balance of MBr in the earth's atmosphere. Based upon

information from several sources, atmospheric release can be broadly attributed to

4 major sources i T a~.e I- 1 with agricultural use of MBr being attributed to 25 % of

the atmospheric MBr. Recently, Yagi et al. (1993) estimated that agricultural use

contributes 80 % of anthropogenic MBr in the atmosphere. MBr is used industrially

as a methylating agent and a degreaser.

More recent investigations emphasize new facets (Table 5-2) of the

uncertainty of various source contributions and sinks for atmospheric MBr. New

modeling studies based upon recent monitoring data suggest that the ocean

absorbs more MBr from the atmosphere than previously thought. Rapid removal by

the ocean shortens the average estimated life-time of atmospheric MBr to just

0.7 year, much less than an earlier estimate of 2 years (Friebele, 1997).

Table 5-2. A recently revised budget for atmospheric MBr (Yvon-Lewis and
Butler, 1997).

Emissions Uptake
MBr Sources (Gg yi) MBr Sinks (Gg y-)
(Gg y') (G y')
Oceans 56 Oceans 77
Fumigation (soils, durables, 46 Hydrolysis 86
perishables and structures) and uv
Gasoline, leaded 15 Soils 43
Biomass Burning 20 Plants ?
Total 137 Total 206

These new results with regard to the role of oceans in the level of

atmospheric MBr (Friebele, 1997) raise a question about the contribution that

anthropogenic MBr plays in destroying the ozone layer. A newly revised MBr budget

(Table 5-2) presented by Yvon-Lewis and Butler (1997) estimates that sinks

(206 Gg y ') greatly exceed sources (137 Gg y'1). These estimates imply that MBr

should not accumulate in the atmosphere. Obviously, high uncertainty prevails

regarding sources and sinks for the accumulation of atmospheric MBr. The buffering

role of oceans implies that reducing agricultural sources will not affect atmospheric

MBr to any great extent, and yet could cause unnecessary economic stress on the

agricultural community. Thus, the economic value and efficacy of MBr fumigation

to control pathogens in agricultural soils warrant renewed research efforts.

Modeling MBr fate and transport during agricultural use (soil fumigation) is

particularly needed to provide improved insight of system dynamics with regard to

atmospheric emissions resulting from MBr fumigation and to provide guidelines for

future research needs.

Table 5-3. Selected chemical and physical properties of MBr (CH3 Br) (Adapted
from Hacherl, 1994).

Selected Property Magnitude
Molecular weight (g mole ') 94.94
Boiling point (oC) 4.6
Solubility in water (g kg'') at 20 C 13.4
Gas density (g L-1 at 20 C) 3.974
Diffusion coefficient (cm2 s-1 in free air at 20 OC) 0.0903
Henry's Law constant at 20 C 0.251
Vapor pressure (mm Hg at 20 OC) 1420

MBr is an odorless hydrocarbon (CH3Br) except at extremely high

concentrations, where it has a sweet odor. It has a low boiling point (4.6 C) and

thus remains in gaseous form at temperatures normally encountered under field

conditions. The important chemical properties of MBr are listed in Table 5-3. Lethal

toxicity of MBr to humans has been noted at 1600 ppm in air for 10-20 hours, or at

7900 ppm for 1.5 hours (USEPA, 1986). The lowest inhalation level found to cause

toxicity in humans is 35 ppm in air (EXTONET, 1997).

Transport Mechanisms for MBr Fumiaant in Soil

Although MBr is extremely volatile, it is commonly injected as a liquid from a

pressurized cylinder. The formulation injected is typically a mixture of 98 % MBr

with 2 % chloropicrin (Hacherl, 1994). Chloropicrin is added to act as a lachrymatory

warning agent. Transport of gaseous MBr in the soil environment is mainly via

diffusion through interconnected air-filled soil pores. Mass flow of gases in the soil

environment has been proposed by several mechanisms, primarily controlled by

changes in temperature and barometric pressure, effects of wind, and forced

movement due to infiltrating water (Baver et al., 1972).

The diffusion coefficient for MBr in free air has been reported to be

325 cm2 hr' at 20 OC (Siebering and Leistra, 1979). The temperature-dependence of

the diffusion coefficient in free air can be defined by use of Eq. [2.15]. Estimation of

the diffusion coefficient of gases in a porous soil system is done by taking into

account the air porosity and tortuosity of the soil system. The tortuosity can be

determined by use of Eq. [3.28].

Maharajh and Walkley (1973) measured MBr diffusion in water for

temperatures ranging from 9.9 to 34.9 C (Fig. 5-1). An almost linear increase in the

diffusion coefficient was observed with temperature over the range of investigation.

o Maharajh and Walkley (1973)
0.07 Linear fit (r = 0.98) o

S0.06 Figure 5-1. Temperature
dependence of MBr
0.05 diffusion in water (Data
0 from Maharajh and
t.0.04 0
0.04 Walldey, 1973).

0.03 0

10 20 40
Temperature (C)3

Fate of I TEr Fi mi q irr .rl til

The fate of MBr is not only dependent upon transport by diffusion in the

gaseous phase, but also upon partitioning between the gaseous and liquid water

phases as well as between the liquid water and solid phases. Reactions in the water

and solid phases greatly influence MBr fate and transport. In general, two ways

exist by which the fate of MBr in soil is affected. The first includes factors that affect

MBr movement in soil in a reversible manner, such as sorption to solids and

dissolution in water. The second includes factors which contribute to the

irreversible chemical degradation or destruction of MBr molecules, such as


Gaseous MBr reversibly partitions into both the solid matrix and the liquid

phase. This partitioning may be described by instantaneous linear sorption under

- ri 1, uir, conditions such as (Brown and Rolston, 1980)

S(R- 1)[5.1]
s = kc ;k, e- p p 5.1
e6 + +p

where k, is a partitioning coefficient between the gaseous phase and the

combination of soil + water phase. Goring (1962) reported the Henry's Law constant

as 0.244 [-] at 20 C for the partitioning of gaseous MBr into liquid water. Mutziger et

al. (1996) reported lab-measured Henry's Law constant values over a temperature

range for the partitioning of MBr into water (Fig. 5-2).

0 Goring (1962)
S0.8 o Observed (1.2 atm)
-- Estimated (1.0 atm) Figure 5-2. Laboratory-
Smeasured and estimated
S0.6 Henry's Law constants at
0 different temperatures for
u 0.4 O the partitioning of MBr
into water (Adapted from
202 Mutziger et al., 1996).

0.0 10 20 30 40 50 60 70
Temperature (oC)

Reversible partitioning reactions are characterized by the fact that, as

gaseous phase MBr concentrations decrease, MBr molecules will desorb from soil

solids or volatilize from water and thus be returned to the gaseous phase.

Irreversible chemical or biological degradation reactions, on the other hand, are

ones in which MBr is converted into another compound in such a way that it cannot

return to the gas phase and thus is permanently lost from the system.

In an irreversible hydrolysis reaction, MBr is known to react with hydroxide

ions (OHE) in aqueous solutions through both unimolecular and bimolecular

nucleophyllic substitutions to form methanol and the bromide (Br) ion as shown in

the following reactions (Moje, 1960)

SNl: H3C:Br- [H3C] +:Br
[H3C] + : OH H C: OH

S2: H3C: Br+ : OH- [HO: CH3: Br] H3C: OH : Br 5.3]

The rate of these hydrolysis reactions is greatly enhanced in the presence of

ultraviolet light, although light from the visible spectrum has no effect (Castro and

Belser, 1981; and Gentile et al., 1989, 1992). Increased pH and increased

temperatures also show an enhancing effect upon the hydrolysis of MBr

(Hacherl, 1994).

Moje (1960) indicated that unimolecular nucleophyllic substitution should be

the major mechanism for the hydrolysis of MBr in water. Thus, the hydrolysis

reaction becomes first-order. Reactions of MBr in soil water with soil organic matter

involve the transfer of methyl groups to the organic matter (Maw and Kempton,

1973). These reactions are expected to be first-order due to the large excess of

organic matter reaction sites (Brown and Rolston, 1980). Brown and Rolston (1980)

measured the appearance of Br in -hl -.i ,, ii ,jiuo rr._- n :i ri nr from the decomposition

of MBr as a function of time for three different soils. The production of Br for a given

introduction of MBr was greatest for organic soil, intermediate for loam soil, and

lowest for sand. The results thus suggest that mechanisms of decomposition

involving clay particles and soil organic matter are significant to the decomposition

of MBr in soils during fumigation. A first-order hydrolysis rate coefficient of

0.014-hr' was reported for Yolo loam soil by Rolston and Glauz (1982). Hornsby et al.

(1996) reported the soil organic carbon sorption coefficient (k_) for MBr to be

22.0 ml g' of organic carbon. In a recent study by Jin and Jury (1995), on a

Pachappa fine sandy loam soil with 0.49 % organic carbon, the degradation rate of

MBr was observed to be first-order and was insensitive to changes in soil water

content (Table 5-4). However the MBr degradation rate showed a slight increase at

water saturation. These results confirmed earlier work showing that organic matter

is the major factor influencing degradation of MBr in soils.

Table 5-4. Measured degradation rate coefficient and half life of MBr as a function
of soil water content (Adapted from Jin and Jury, 1995).

Degree of Water Saturation in % Half-life Rate coefficient
(Water Content, 0,) (days) (hr-')
45 (6, = 11.5) 12.5 0.0023
76 (9L = 19.4) 12.5 0.0023
100 (6, = 25.6) 10.7 0.0027

Although hydrolysis deactivates toxic MBr molecules, toxic levels of bromide

may accumulate in the soil water of treated soils (Macartney and Price, 1988; and

Maw and Kempton, 1982). High levels of bromide accumulation have also been

found in crops such as tomato, carrot, lima bean, lettuce, radish, tobacco, wheat,

orange, and lemon grown on initially fumigated soils (Hoffman and Malkomes, 1979;

and Maw and Kempton, 1973).

Permeation/Transmission of MBr through Polyethylene Mulch

Open-field fumigation without a plastic soil cover (mulch) may provide a

large input of MBr into the atmosphere. In plastic-mulch culture, atmospheric

emission occurs by diffusive transport (or permeation) through the plastic mulch

(within the bed) as well as by diffusion through uncovered soil surfaces in furrows

between parallel soil beds. Several environmental factors contribute to diffusive flux

through the soil in furrows and through the plastic mulch in beds. Environmental

factors with the greatest impact upon diffusive flux include the temperature

gradient across the plastic mulch, the moisture or condensed water layer on the

underside of the plastic mulch, and the sweeping effect of wind over the surface

(overland wind velocity).

The physical dynamics of water and heat transport beneath plastic-mulched

soil surfaces are quite complex (as discussed in Chapter-2.). Temperature increases

at the plastic mulch during the daytime by as much as 15 to 20 oC due to solar

radiation are not uncommon (Ham and Kluitenberg, 1994). Such temperature

increases during daylight hours enhance the water evaporation rate, thereby

increasing relative humidity of the air gap and of soil beneath the plastic mulch.

Since the plastic cools faster than the soil, the condensation of water from the

humid soil air tends to occur on the underside of the plastic mulch during nighttime

hours. As a result, permeability of the plastic mulch to MBr would be reduced due to

a blocking effect by the layer of condensed water. As the mulch temperature

increases during daytime hours, evaporation then removes the condensed water

layer. In cases where permeating MBr is soluble in or sorbs to the plastic material,

the rate of movement by diffusion would further be reduced, thus delaying the

equilibrium. As such, only operational value for the diffusion coefficient can be


Table 5-5. Ratet of escape of MBr through several films at various temperatures
(Adapted from Kolbezen and Abu-El-Haj, 1977).

Temp. LDPE* HDPE-l* HDPE-2" Polybutene Saran Black LDPE
(C) (1-miil) (-mil) (1-mil) (1-mil) (0.5-mil) (6-mil")
23 8.2 2.7 1.4 2.9 Nil 1.6
30 11.1 3.5 1.7 4.0 Nil 2.2
40 15.1 5.0 2.9 6.3 Nil 3.4
50 20.1 6.8 3.9 9.1 0.18 5.9
60 30.6 8.8 5.2 11.8 0.41 9.2

t Units are ml hr' -2 m-2 L',
low density polyethylene,
" high density polyethylene sample 1,
tt high density polyethylene sample 2, and
* 1 mil =2.544E-3 cm.

Few investigations have been conducted to measure the diffusion coefficient

of MBr through polyethylene sheets used for mulching. Kolbezen and Abu-El-Haj

(1977) measured rates of escape of gaseous MBr through six selected films at

different temperatures. The escape rates increased with temperature for all films.

MBr escape rates ranged from 0.41 to 30.6 ml hr-' m-2 ml' L-' for the six materials.

The authors found one of the high-density polyethylenes (Table 5-5) to be an

effective barrier to MBr transport, with minimum escape rate over the temperature

range investigated.

During a field experiment for unbedded mulched soil in California Yagi et al.

(1993) observed that the cumulative emission of MBr decreased as the thickness of

LDPE plastic increased (1.03,1.29,1.37, and 1.6 mil thickness were used).

Unpublished results from field investigations for mulched soil bed system at the

University of Florida Green Acres Experimental Station, near Gainesville, FL, with

black LDPE (1.2 mil) and clear laminated plastic (relatively thicker than black

LDPE) showed emission fluxes differed between the two plastic materials.

Emissions were reported for both the plastic-covered beds and uncovered furrows in

between the parallel beds.

0.25 Clear Laminated Plastic
0.20 Feb. 1996 Bed: Observed
0.15 Bed: Moving Average
0.10 Furrow: Observed
0.05 "---------- Furrow: Moving Average


0.25 Black Plastic
S( Feb. 1996 Bed: Observed
0.15 \ Bed: Moving Average
n0.10- Furrow: Observed
S ----------. Furrow: Moving Average

a 3.0
" 3 Black Plastic
, June 1995 Furrow: Observed
W 2.0 Furrnw: Movino AvrFur

Bed: Observed
--- Bed: Moving Average

- ..-

Figure 5-3. MBr
emission fluxes in
plastic mulched beds
during fumigation at
the University of
Florida Green Acres
Experimental Station
located near
Gainesville, FL.

24 48 72 96 120 144 168
Time after fumigation (hour)

MBr was applied at the rate of 450 kg ha-' with a 3-nozzle applicator at a

depth of 20 and 33 cm during the summer of 1995 and winter of 1996, respectively,

in soil beds of Arredondo Fine Sand soil used for vegetable production. As expected,

higher magnitudes of the emission fluxes were observed in the summer shallow


injection over winter deeper injection during fumigation (Fig. 5-3). A diurnal cyclical

pattern was evident in the emission fluxes to the atmosphere due to non-isothermal

conditions. The relatively thicker clear laminated plastic was able to retain slightly

more MBr in the soil bed over black LDPE (Fig. 5-3). However, material balance

estimations from emission fluxes and soil Br residue measurements for MBr applied

during fumigation revealed that 39, 46, and 70 % of applied MBr were

unaccountable in bed systems with black LDPE (June, 1995), black LDPE (Feb.,

1996), and clear laminated plastic (Feb., 1996) mulches, respectively. Hence, definite

conclusions were not feasible from these field investigations of MBr fumigation in

plastic mulched vegetable beds (personal communication with investigators').

Hacherl (1994) conducted laboratory investigations measuring the diffusion

through black LDPE (1.25 mil) film commonly used in agriculture. Figure 5-4 shows

the results of diffusion investigations by Hacherl (1994) and Kolbezen and

Abu-El-Haj (1977) for black LDPE with 1.25-mil and 1.00-mil thicknesses,

respectively. The results revealed a near-linear increase of diffusion coefficient with

temperature, with lower values in the case of Hacherl (1994) corresponding to

thicker LDPE. In another study (Jin and Jury, 1995), from the data presented, the

diffusion coefficient was estimated to be 0.005969 cm2 hr- for black LDPE (1.0 mil).

SA. G. Hornsby, L.T. Ou, D. Dickson, S.J. Locascio, and P.S.C. Rao

S0.o8 Kolbezen and Abu-El-Haj (1977)
Black LDPE: 1.00- mil
0.007 r =0.97
0.006 Figure 5-4. MBr diffusion
0 through black LDPE of
S0.005/ 7 different thicknesses (Data
S0.004 from Hacherl, 1994; and
C. Kolbezan and Abu-El-Haj,
r 0.003
003 o 1977).
0.002 e Hacherl (1994)
So-- 0 Black LDPE: 1.25-mil
0.001 r = 0.99
20 30 40 50 60
Temperature (C)

Modeling of MBr Transport during Fumigation

Only a few prior modeling efforts have been attempted to describe MBr

transport in soil during fumigation. Hemwall (1959) described a two-dimensional

model for gaseous diffusion of MBr in soil with relatively low water content. The

model assumed a small constant water content and isothermal soil conditions.

Typical simulations were presented to demonstrate time-dependent patterns of

fumigant concentration and biological control within an uncovered soil profile after


Siebering and Leistra (1979) presented a one-dimensional model for gaseous

diffusion of MBr under isothermal conditions. Water content was held constant and

small, such that water flow could be ignored. Comparison of model simulations with

experimental results for MBr concentrations during soil fumigation applied to soil

under cover in a greenhouse matched only in the upper part of the soil profile

(to 40 cm) but deviated substantially elsewhere.

Rolston and Glauz (1982) proposed a model describing radial diffusive

transport of gaseous MBr away from injection chisels during fumigation of a field

profile. The model was based on similar work done by Hemwall (1959, 1960) and

Siebering and Leistra (1979). The model included the dissolution-distillation of MBr

gas on soil solids and in soil water, using both kinetic and equilibrium relationships

and hydrolysis to yield Bf ions according to first-order kinetics. The model

considered the case where the soil surface was covered with a barrier totally

impermeable to the diffusion of gas along with the case of a bare soil surface.

Reported simulation results provided a better match with the experimental data of

Abdalla et al. (1974), which involved a mulch barrier, when no barrier was

considered at the soil surface during simulations. Sensitivity analysis of input

parameters showed reaction kinetics to be unimportant, so they could be safely

ignored. Water flow and non-isothermal temperature effects were ignored in all

simulations. Also, in that model, the boundary conditions (BC) at the soil surface

were unrealistically handled by making the plastic mulch impermeable and

applying a zero-concentration fixed BC at the soil surface.

Mignard and Benet (1989) considered a kinetic relationship for MBr

partitioning between the gaseous and solution phases. Synthetic simulations were

conducted and compared with equilibrium models, and profiles were developed to

forecast the treatment zones and behavior of MBr. The authors inferred from their

analysis that a kinetic partitioning relationship is important only for rapidly

changing boundary conditions. This model also ignored water and temperature

variations in the soil.


Recently Wang et aL (1997) adopted a two-dimensional (2-D) finite-element

model, CHAIN_2D (Sim~mek and van Genuchten, 1994) to simulate the fate and

transport of MBr fumigant in soil. CHAIN_2D numerically solves the partial

differential equations for 2-D nonlinear, nonequilibrium simultaneous transport of

water and heat, (but not coupled heat-water) and solute in a variably saturated

porous medium. The temperature dependence of various coefficients (such as DI,

K,, etc.) were described in similar fashion to the present work. However, the model

lacked a realistic description of the soil-surface plastic-mulch/atmospheric

boundary. Placement of a plastic mulch significantly changes the energy balance at

the soil/atmospheric interface (as discussed in Chapter-3) and requires proper

redefinition of the energy exchange and balance at the interface to simulate non-

isothermal soil conditions. The transfer of MBr to the atmosphere was described

(Wang et al., 1997) by introducing a transfer coefficient (h,) as

FlCux =--- (g-g)= hr, (g-g.) [5.41

where db (cm) was considered a composite boundary thickness inclusive of an air

gap below the mulch, the mulch thickness, and an overlying atmospheric boundary

layer; and where h, was described as a transfer coefficient (cm s '). The db

parameter was shown to be temperature-dependent through use of an equation

similar to Eq. [2.14], but the scientific basis for doing so was not explained.

Although a cyclic nature in emission flux in accordance with non-isothermal

conditions was demonstrated, model simulations failed to accurately describe

(Fig. 5-5) the cyclic field-measured emission flux reported by Yates et al. (1996).

120.0 9 -- Observed
100.0 -


40.0 -

20.0 0 0
0.0 0 6 0 .
.0 10 20 30 40 50 60 70 80 90 100 110 120 130
Time (hour)
Figure 5-5. Simulated and measured flux density through plastic-
mulched surface. Adapted from Wang et al. (1997).

It is thus evident that, in modeling MBr fate and transport, various

investigators have made simplifying assumptions (such as small and constant soil

water content and isothermal conditions) that reflect unrealistic field conditions.

Also, BCs at the soil surface which are more realistically described by a Neumann-

type in the case of MBr fumigation have been improperly designated (such as a zero

concentration fixed BC at nodes exposed to the atmosphere for uncovered soil). The

dissolution (Henry's Law constant) and diffusion (both in air and water) of MBr are

temperature-dependent, and variably saturated soil-water flow is a common

occurrence under field conditions. This is particularly true for humid climates.

Analysis of Fate and Transport of MBr Fumicant

The model presented in this dissertation represents many improvements over

earlier efforts (Hemwall, 1959; Mignard and Benet, 1989; Rolston and Glauz, 1982;

and Siebering and Leistra, 1979), for it includes coupled transient flows of both heat

and water based upon PDT theory. Realistic plastic-mulch boundary conditions are

described thoroughly by inclusion of optical properties of the plastic mulch in the

energy balance at the soil/atmospheric interface. This model provides opportunity to

consider the effect of improved temperature upon estimates of solute transport (both

in the liquid and gaseous phases).

Modeling investigations (synthetic simulations) of the fate and transport of

MBr fumigant were conducted to obtain improved insight into the system dynamics

for MBr fumigation of soils. Modeling simulations of MBr fumigation prior to this

research have utilized overly simplified assumptions, such as gaseous diffusive

transport in air-filled soil pores under isothermal conditions with the effects of

transient liquid-water flow neglected. Simulations reported here utilized an

assumption that was more realistic of actual field conditions; i.e., non-isothermal

soil conditions that fluctuate diurnally due to net solar radiation, thus altering vapor

and liquid water flow as well as MBr diffusion. Objectives were as follows:

1) To investigate important major processes influencing MBr fate and transport

including volatilization and degradation in fumigated soil beds used for

vegetable and fruit production in Florida,

2) To investigate the effect of non-isothermal and variably saturated soil

conditions on material balance of MBr during simulated fumigations,

3) To evaluate the effectiveness of the current management technique-i.e.,

plastic-mulching of soil beds-in minimizing MBr emission from fumigated

soil beds to the atmosphere, and

4) To examine alternative management methods for decreasing MBr emission

during commercial use of soil fumigation.


Simulation analysis was designed considering major factors known to affect

the fate and transport of MBr. The following factors were considered in the design of


I. Effect of variable water saturation in the soil bed: Two levels of heterogenous

water saturation were selected to represent a relatively dry soil and a

reasonably wet profile. The effects of hydrolysis (degradation) on MBr

material balance were considered for different levels of soil water saturation.

II. Effect of variable temperature regime in the soil bed:

A. Non-isothermal conditions: Effects of transient thermal conditions

were considered upon the following input parameters for the model

1. MBr diffusion coefficient in the air phase (DI),

2. MBr diffusion coefficient in the water phase (D"),

3. MBr partition coefficient between gaseous and aqueous

phases (Henry's Law constant, K,),


4. Diffusion coefficient of MBr through polyethylene film (D"fpc),

B. Isothermal conditions: Uniform thermal conditions were provided for

comparison (average daily temperature was assumed to exist

uniformly over the simulation domain).

III. MBr diffusion through plastic mulch: Due to uncertainty in published data

and based on reported field investigations of MBr emissions from plastic-

mulched soil, different magnitudes for the MBr diffusion coefficient through

plastic film were used to examine the material balance of MBr.

IV. Depth of fumigant injection: Two depths of fumigant injection zones were


A. Typical shallow injection zone: 33 cm depth total

B. Deep injection zone: 66 cm depth total (arbitrarily chosen for


V. Effect of irrigation (dousing of water on the bed): Effect of irrigation was

examined to evaluate the surface-sealing effect (due to wetness) on

reduction of gaseous emissions for MBr.

VI. Fumigation Efficiency: Soil sterilization during fumigation was investigated

based on the nematoxic level of MBr (20 ppm or greater) and exposure time

(minimum 24 hr) for 2nd-stage juveniles of root knot nematode (Mckenry and

Hesse, 1978).

VII. Hydrolysis of MBr: As no data for the hydrolysis rate of MBr in Florida soils

was available, sensitivity analysis was carried out for 3 rates; 1) 0.014 hr' (for

California, Yolo loam soil, Rolston and Glauz (1980)), 2) 0.007 hr ', and

3) 0.0035 hr '.

VIII. Diurnal water and heat dynamics in soil bed: The influence of atmospheric

interactions with the soil surface (both bare furrow and mulched bed) were

examined by analyzing the water and heat dynamics in the soil bed system

over a diurnal period of 24 hrs.


Arredondo Fine Sand (Grossarenic Paleudults, loamy, siliceous,

hyperthermic), a coarse-textured soil at the IFAS Green Acres farm near Gainesville

and commonly found in north central Florida, was selected for the modeling

investigation. Representative soil properties were taken from 'Soil Science Research

Report # 88-1, 19881. The soil consisted of five horizons namely; Ap (0-23 cm),

El (23-66 cm), E2 (66-99 cm), Btl (99-165 cm), and Bt2 (165-203 cm). The weather

data for simulation purposes were taken from a weather station at the Green Acres

Farm of IFAS located west of Gainesville. The month of February was selected for

all simulations, since MBr fumigations are normally performed for the Gainesville

area at that time of the year.

A spatial profile discretization of the cross-sectional soil bed resulted in 2782

elements with 1474 nodes (Fig. 5-6). Bed symmetry was assumed, so only a half-

section (AB-BC-CD-DE-FE-AF) of the soil bed needed be considered for simulations.

'Compiled by Univ. of Florida, IFAS, Soil Science Dep. Soil Characterization Lab. in
cooperation with the USDA-SCS.

No flow BCs on the left- and right-hand sides (AF and DE) were defined for heat,

water, and solutes. The DE boundary represents a vertical symmetry line through

the middle of a furrow between two parallel soil beds. The AF boundary represents a

vertical symmetry line through the middle of a soil bed. The lower boundary (FE)

was designated as a zero thermal gradient with gravitational water flow (= K(h, T)).

-x Mulched-Bed

-. -- ---- .. .. ..
Dis c fr:- r Be
i : ... i .:_ ~ :3 .:
.Inection Zre. I .

-75- : 2

-. 7-- ------- ---.
B l .

4I ... ... ....- .

-175 -
t B3 2 .

I[I l I: .:. I.I

Distance from Bed Center (cm)

Figure 5-6. Schematic of the discretized simulation domain of
the soil bed.

For solute transport, both convective (in the liquid phase) and diffusive (both

in the aqueous liquid and gaseous phases) transport were allowed at the bottom

(FE), with the assumption that solute leaving the lower boundary is lost from the

flow domain (i.e., there is no return flow). The emission loss of MBr to the

atmosphere (AB-BC) was described through diffusive flux, assuming a 1-cm

boundary layer and zero concentration on the atmospheric side of the boundary

layer. These emission fluxes may be defined explicitly as

a) For Bed

'd = gasi s gs+ l(air)
S I Zs Zs+ (air) /

b) For Furrow

jFwrow i 8 s(air) +,ni g -l
Jgas free air --, +ZI- O),C z--- z- j

where the subscripts s, +,, and, represent nodes at the soil surface, above the

surface, and below the surface, respectively.

Table 5-6. Physical properties for 5 distinct horizons of a profile of Arredondo fine
sand soil.

Soil Layer Sand Silt Clay OCt BD K at
Soil Layer (%) (%) (%) (%) (g cm 3) (cm hr
Ap 91.9 3.3 4.8 0.83 1.52 6.2
El 92.0 3.3 4.7 0.28 1.54 16.7
E2 91.4 3.0 5.6 0.14 1.52 15.1
Btl 86.0 3.7 10.3 0.11 1.56 6.4
Bt2 75.2 3.0 21.8 0.20 1.57 2.4

SOC- Organic Carbon,
BD- Bulk Density

The data for physical properties of the soil profile used are presented in

Table 5-6, with other related parameter values in Table 5-7. The soil water

characteristic data used for estimating (van Genuchten, 1980) hydraulic-function

parameters were taken from the report mentioned earlier. Average weather

parameters' used in simulations for the month of February were: 1) average daily

temperature 11.6 C, 2) temperature amplitude 7.7 OC, 3) wind velocity 11x104

cm hr1, 4) day-length 11.0 hrs, and 5) global daily radiation 1400 J cm-2 day'.

Table 5-7. Parameter values for 5 distinct horizons of a profile of Anedondo fine
sand soil.

Soil pI I c K [ 1 Kd2 van Genuchten Parameters
Layer (hr-1) (cm) (cmgg-) (cm ( g c) a (cm 0
Ap 0.014 0.1826 0 0.0225 2.124 0.43 0.03
El 0.014 0.0616 0 0.0216 2.762 0.42 0.04
E2 0.014 0.0308 0 0.0271 2.626 0.43 0.04
Btl 0.014 0.0242 0 0.0211 2.091 0.41 0.06
Bt2 0.014 0.0440 0 0.0307 1.557 0.41 0.13

t- hydrolysis rate coefficient for MBr,
soil sorption coefficient for MBr, and
tt- soil sorption coefficient for bromide ion.

Temperature-dependence of the diffusion coefficients for MBr in gaseous

and aqueous phases of the soil were defined by use of Eq. [2.15] and data reported

ir, F-i. 5-1 from Maharajh and Walkley (1973), respectively. As information is

unavailable for the temperature dependence of the hydrolysis rate coefficient lp,), a

fixed value of 0.014 hr-1 was used (Rolston and Glauz, 1982). The diffusion coefficient

'Climatological Data, Agronomy Department and NOAA Cooperating, 1996.

for the plastic mulch was defined using the data of Hacherl (1994). The MBr

fumigant was assumed to be applied at a rate of 450 kg ha-' with a three-shank

applicator, so that 1/ shank (3591 pmoles) injections occur within the flow domain

shown in Fig. 5-6. Unless otherwise noted, the depth of the injection and time of day

for injection for all the simulations were set at 33 cm and 9:00 AM, respectively.

Black LDPE of thickness 1.25 mil was assumed as the plastic mulch m all the

simulations. Mass and energy balance errors were less than 1.0 % for all simulations

reported in this work.

Results and Discussion

Plastic-Mulching of Soil Beds to Control MBr Loss to the Atmosphere

Preliminary simulations conducted using the laboratory-measured data

(Fig. 5-7, Level-1) of Hacherl (1994) for diffusion through plastic (black LDPE,

1.25-mil thick) provided estimates for the atmospheric losses of MBr during

fumigation of the soil bed that were very low (23 % during the 1" day) compared to

reported values of 40 to 60 % for recent investigations (Yagi et al., 1993; Yates et al.

1996). Hence, parallel elevated temperature-dependent levels were designed

(Fig. 5-7) for investigation. The use of an experimentally measured diffusion

coefficient for a 1.25-mil thick black LDPE by Hacherl (1994) in simulations showed

that very little MBr escaped (0.02 % of that applied) through the plastic mulch into

the atmosphere (i.e. from the bed excluding the furrow) during the 1st day. The

diffusive flux (Fig. 5-8: Level-1, and Fig. 5-9: a, Level-1) based upon the measured

diffusion coefficient of plastic mulch from the soil bed was very low and it was also

27 Elevated Level-4 (= Level-I x 104)
22 6430*s^
12 _
S 2.7 Elevated Level-3 (= Level-I x 103) A

1' 2.2 1.2 -._

8 0.27 Elevated Level-2 (= Level-1 x 102)
S0.12 _______________
0.0027 Measured Level-I (Hacherl, 1994)
0.0012 20
20 25 30 35 40
Temperature (C)

Figure 5-7. Measured and elevated levels of diffusion coefficients
through 1.25-mil black LDPE plastic mulch used in the
model simulations.

demonstrated by the material balance after 1 and 3 days (Table 5-8) that only 23.4

and 45.8 % of the applied MBr, respectively, was lost from the system (mulched soil-

bed plus bare furrow) to the atmosphere. This result is much lower than the 40 to

60 % atmospheric losses earlier reported from field investigations (Yagi et al., 1993;

Yates et al., 1996), which were carried out under unbedded conditions with

complete mulch cover of the soil. Initial average fluxes both from the bed and the

furrow were substantially higher than for later periods (after 24 hours) when the


concentrations in soil had decreased due to spreading, hydrolysis, and loss of MBr to

the atmosphere. Also note that emission fluxes from the furrow were nearly 5 times

larger than for the bed, except for elevated plastic-mulch diffusion Level-4 due to

the absence of a mulch barrier over the furrow.



35 00
8 0 7_ ----------------


\ 5.00

I 3A0 2


2, 125 ----- Elevated Level-2
S-- Elevated Level-3
1.50 -- Elevated Level-4
Time after fumigation (hour)

Figure 5-8. Losses of MBr to the atmosphere through soil bed and furrow
surfaces using measured and elevated diffusion coefficients for
MBr through plastic mulch'.

'Small steps in some of the figures resulted from the need for plotting more points
within certain time periods.


Use of elevated levels of diffusion coefficients for plastic mulch showed more

losses of MBr into the atmosphere from the system (Level-2 = 25.1 %; Level-3 =

34.7 %; and Level-4 = 51.7 %) one day after fumigant injection. In all cases,

emission fluxes from the bare soil furrow were highest near the bed due to close

proximity to the injection zone. Figure 5-9: a shows MBr fluxes in the bed 6 hours

after fumigant injection, when peak emission fluxes occurred from the bed (Fig. 5-8).

(a) Measured, Hacherl (1994) Elevated Level-2
.. .I, m .I .

Elevated Level-3 Elevated Level-4
t '- 1 .--.

-- '*2-- .*0

S 15 30 45 60 75 0 15 30 45 60 75
Distance from bed center (cm)

Figure 5-9. MBr total concentration (pmol cm- ) contours and gaseous
MBr flux vectors (pmol cm2 hr-') in the bed 6 hrs after infection
of MBr fumigant.
a) 4 levels of diffusion coefficients through the plastic mulch,
b) non-mulched bed.

Emission fluxes (Fig. 5-9: a) increased from the mulched bed with elevated

levels of diffusion coefficients for the plastic mulch. Also to be noted here is the fact

that the trend of diffusive flux was reversed in the furrow compared to that for the

bed, and lower diffusion coefficients for plastic mulch showed higher diffusive flux

and more cumulative loss from the furrow (Fig. 5-8 Fig. 5-9: a, and Table 5-8).

Figure 5-9: b shows the emission of MBr in furrow and bed without a plastic mulch,

for comparison purposes. The emission fluxes were higher from soil-surface regions

close to the fumigant injection zone (Fig. 5-9: b).

(b) No Mulch
0 ltt s-l l:1 lmol cm2 hr-'

Fiqure 5-9--continued

-75 '

-0 15 30 45 60 75
Distance from bed center (cm)

Overall, total material losses of MBr from the simulated soil-system to the

atmosphere (Table 5-8) increased when elevated levels of diffusion coefficients for

plastic mulch were used in the simulations. Three-day simulations were

appropriate, as the results showed that most of the MBr losses occur within the

initial 24 to 30 hrs after fumigation (Fig. 5-8 and Table 5-8). Atmospheric losses were

observed to continue in small amounts even after 3 days, as long as undegraded

MBr is present in the soil.

Obviously much uncertainty prevails regarding the diffusion coefficient of

MBr in plastic mulches. The use of a lab-measured diffusion coefficient for

simulating field conditions greatly under-estimated published results from the field

investigations done in the past. However, the use of elevated diffusion coefficients

(Level-2, Level-3, and Level-4) improved simulation results. It was, therefore,

decided to use the Level-3 diffusion coefficient for all other simulations, as results

obtained from this level were close to those from the published field investigations.

Table 5-8. Material balance (unit: %) of MBr under measured and enhanced levels
of diffusion coefficient for plastic mulch.

Atmospheric Emission
Day-1 Total MBr Losses Hydrolyzed Deep Down-
in Profile flow
Bed Furrow Total
Level-1 69.55 0.02 23.39 23.41 7.04 0
Level-2 67.91 1.81 1 23.29 25.10 6.98 0
Level-3 58.54 13.15 21.57 34.72 6.73 0
Level-4 42.24 34.24 17.41 51.65 6.11 0
No-Mulch 37.29 56.92 5.79 0
Level-1 36.38 0.20 i 45.63 45.83 17.79 0
Level-2 33.94 4.24 I 44.42 48.66 17.39 0
Level-3 23.46 24.55 1 36.52 1 61.07 15.46 0
Level-4 15.05 46.32 1 25.84 1 72.16 12.79 0
N h 13.40 I74.68 111
No-Mulch 13.40 1 74.68 11.91 0

Non-Isothermal versus Isothermal Soil-Conditions

An isothermal simulation was performed in order to compare with

non-isothermal simulations existing under a black LDPE. Thermal effects on solute

transport (involving D9, D",kg, Dgp ,c parameters) were ignored in isothermal (with

respect to solute transport only) simulations. However, the coupled heat and water

flow model was used in its original form, with solute transport being made

independent of temperature fluctuations. The parameters for solute transport were

estimated and kept constant for a daily average temperature of 11.6 C, -.p,,: ji ,

existing during the month of February in the vicinity of Gainesville (diffusion

coefficient in air = 308 cm2 hr-1, Henry's Law constant = 0.15 [-]).

S5.00 Non-Isothermal
S3.75 ----- Isothermal

o Furrow

R 0.99


80.33 Bed
CO~ ~ ------------^
< "000 12 24 36 48 60 72
Time after fumigation (hour)

Figure 5-10. Emission fluxes of MBr through bed and furrow under
isothermal and non-isothermal transport of fumigant.


The isothermal and non-isothermal simulations conducted for a 3-day period

showed (Fig. 5-10) that MBr emission fluxes from both bed and furrow were affected

by non-isothermal conditions. Figure 5-10 shows the diurnal fluctuations in MBr

emissions from bed and furrow in the case of non-isothermal simulations for black-

LDPE, whereas a similar cyclical nature of MBr emission was missing m isothermal

simulations. This emission pattern has been reported earlier (Fig. 5-5) for field

investigations (Yagi et al., 1993; and Yates et al., 1996). Higher diffusive flux rates

for MBr were observed initially from the furrow, which subsequently flattened out

due to lower concentration gradients as increasing amounts of MBr were lost to the

atmosphere. The delay in reaching the initial peak emission flux (Fig. 5-10) from the

bed occurred due to a build-up time for MBr concentrations beneath the mulch and

an increase in temperature of the plastic mulch as the day progressed (application

time 9:00 AM), whereas the bare furrow surface offered a quick passageway to the


Table 5-9. Material balance (unit: %) of MBr for fumigation under isothermal and
non-isothermal conditions.

13.44 21.59 34.73

Total MBr Atmospheric Emission Losses
in Profile Bed Furrow Total

Isothermal 61.28 11.42 20.33 31.75
i j




The results suggest that inclusion of non-isothermal conditions is essential

to correctly estimate the fate and transport of volatile chemicals. Non-inclusion of

temperature effects on solute transport may underestimate (Table 5-9) matenal

losses to the atmosphere in the case of volatile chemicals such as MBr. The

isothermal simulation estimated 260.2 kg ha'1 (57.8 % of application) emission loss

to atmosphere compared to 274.8 kg ha' (61.1 % of application) as a non-isothermal

estimate over a 3-day period.

Variable Water Saturation of Soil

Volumetric water contents present in the soil bed at the time of fumigation

not only affect the hydrolysis of MBr but also its transport in the soil. Two different

wetting regimes were selected for the simulations. An initially saturated soil profile

was allowed to drain and evaporate water for 40 days. Two wetting levels were then

selected to provide the soil water regime:

1) Level-1 (WL-1): 3292 cm3 total water content in the flow domain, and

2) Level-2 (WL-2): 1864 cm3 total water content in the flow domain.

The simulations revealed (Fig. 5-11) that larger amounts of MBr were lost to

the atmosphere when there was less water (WL-2
result was expected, due to less hydrolysis and enhanced diffusion. Higher

temperatures and lower water contents existed (Fig. 5-12) near the surface in the

case of the drier soil system (WL-2). Greater diffusion pathways (larger diffusion

volume due to lower water contents) at higher diffusion rates also existed in the

drier profile. This resulted in a higher MBr emission flux from both the bed and


furrow in the case of the drier (WL-2) profile (Fig. 5-12). The emission fluxes from the

bed were also higher (Fig. 5-11) during the later periods (after approximately 34 hrs)

in WL-1. This was due to the higher water contents in WL-1 causing larger amounts

of MBr to partition into the aqueous phase, which were then released to the gaseous

phase as atmospheric losses occurred in order to maintain equilibrium between the

aqueous and gaseous phases (controlled by the Henry's Law constant).

5.00 WL-1
S3.75 ---- WL-2

S2 50 Furrow

o o.oo '1/-----------
1 .25 1- --------------------------
0.99 -

Time after fumigation(hour)
'S 0.66 -
S ;Bed


< .00^ -^-^ --^ -- -- -*f
Time after fumigation(hour)
Figure 5-11. Emission fluxes of MBr through bed and furrow dunng
fumigation performed under different profile
water-saturation scenarios.

Figure 5-12 shows the plots of water content, temperature, MBr

concentration, and MBr fluxes in the bed under different water saturations (WL-1

and WL-2) for 0.5, 3, 6, 12, 18, and 24 hours after fumigant injection. Temperatures

under the bed and at the bare soil surface were higher in WL-2 as compared to

WL-1, whereas water contents were lower in WL-2 as compared to WL-1 at the

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EW1EGXD3U_T33EQT INGEST_TIME 2013-01-23T14:21:43Z PACKAGE AA00012978_00001