Water and nutrient movement in two tropical cropping systems


Material Information

Water and nutrient movement in two tropical cropping systems
Physical Description:
vi, 170 leaves : ill. ; 28 cm.
Seyfried, Mark S., 1954-
Publication Date:


Subjects / Keywords:
Soils -- Costa Rica   ( lcsh )
Cropping systems -- Costa Rica   ( lcsh )
Soil Science thesis Ph. D
Dissertations, Academic -- Soil Science -- UF
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1986.
Bibliography: leaves 157-168.
Statement of Responsibility:
by Mark S. Seyfried.
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 - 000879791
oclc - 15096901
notis - AEH7576
sobekcm - AA00004870_00001
System ID:

Full Text








In working on this project I have been extremely

fortunate to have had a great deal of assistance. At the

outset, the grant money and connections supplied by Bob Volk

enabled me to get my "foot in the door" in Costa Rica.

While I was in Costa Rica, Carlos Burgos of CATIE looked

after my interests in every conceivable way. He not only

contributed experimental plots, but helped me obtain fund-

ing, allowed me access to transportation and provided the

services of his assistant, Carlos Arya, who was very help-


I owe a debt of gratitude to several others at CATIE.

Roberto Diaz gave me the run of his laboratory and saw to it

that his assistants, Taco, Flaco, and Eduardo, were able to

help me when it was needed. Louis Alpizar provided help

with background data. Gustavo Enriquez provided land for

plot work, and Raul Moreno helped me obtain employment


While in Costa Rica I also received help from a number

of Floridians. In particular, Bob Mansell saved the project

by supplying a much needed neutron probe. Jack Ewel re-

lieved me from the daily routine long enough to take a short

vacation and shared data with me. And Chris McVoy assisted

with interest, discussion and companionship.

In Florida, Suresh Rao "adopted" me as his student,

which is probably the best thing that could have happened.

His input into all phases of the research, even in his

absence, has been invaluable. Ron Jessup has also taught

and helped a great deal. I am especially grateful to these

two men.

My remaining committee members, Don Graetz, Nick

Comerford, and Jerry Bennet have provided valuable feedback

and definitely improved the end product. In addition, I

must thank Peter Nkedi-Kizza and Linda Lee.

Finally, last and most, I thank my wife, Helen Fisher,

who not only put up with a great deal, but helped me most

when I most needed help.





ABSTRACT..................................................... v



Leaching in Humid Tropical Regions..................3
Presentation........................................ 8


Introduction......................................... 9
Materials and Methods..............................14
Results and Discussion.............................. 21


Materials and Methods..............................62
Summary ...........................................105


Introduction.......................... ............108
Materials and Methods.............................115
Discussion......................................... 133
Conclusions....................................... 143


Water Movement and Nutrient Leaching...............146
Mineralization and Nutrient Cycling................150
Overall Assessment of Cropping Systems.............152
Future Work........................................154


BIOGRAPHICAL SKETCH......................................169

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



Mark S. Seyfried

August 1986

Chairman: P.S.C. Rao
Major Department: Soil Science Department

Soil-water and nutrient movement in two contrasting

cropping systems was compared to investigate the effect of

cropping system management on leaching losses under humid

tropical conditions. One cropping system, the monocropped

annual, consisted of maize (Zea mays L.) alone; the other,

mixed perennial system, consisted of cacao (Theobroma

cacao), laurel (Cordia alliodora), and plantain (Musa

paradisiaca) grown together. A model was developed to

calculate nutrient leaching losses given measured soil

solution concentrations. Displacement through undisturbed

soil columns was used to select a modeling approach. The

contribution of mineralized soil N was studied in incubation


Miscible displacement studies indicated that nutrient

movement was well described by the convective-dispersive

model except when the soil was at or very near saturation.

Saturated hydraulic conductivity values measured in the

field were considerably higher than rainfall intensities,

suggesting that preferential flow conditions were very rare.

A zero-order kinetic model, with consideration made for

pretreatment effects, was shown to describe the measured N

mineralization rates well. These results appear to apply to

a variety of soils.

The simulation model used was based on soil-water

capacity parameters. Piston displacement of solutes and

drainage to field capacity were assumed.

Simulated profile soil-water contents closely approxi-

mated those measured. Calculated net leaching losses of N,

K Ca2+ and Mg2+ were greater in the monocropped annual

system due mainly to differences in soil solution concentra-

tions, which were constant over time. The calculated loss

of N from the monocropped plot was 56 kg ha-1 over 260 days.

Losses from the mixed perennial plot were negligible.

The solute residence time within the crop root zone was

proposed as an index of cropping system sensitivity to

leaching loss. The residence time in the mixed perennial

system was about 1.4 times larger in the monocropped system.

The low leaching losses from the mixed perennial plot, in

spite of substantial annual inputs of nutrients, indicates

that that system is conservative in terms of nutrient use.

This was attributed to the relatively long residence time in

the rooting zone of the mixed perennial plot.


The humid tropics, as defined by Sanchez et al. (1982),

comprises about 10% of the earth's land surface. This

region is currently undergoing very rapid ecologic, economic

and social changes. These changes frequently require more

intensive use of the land to produce food and fiber. This,

in turn, requires changes in extant farming methods.

The traditional system of farming for much of the

region, which is still used extensively, is commonly re-

ferred to as slash and burn agriculture or shifting cultiva-

tion. Although the details of operation are extremely

variable, the basic pattern is fairly consistent throughout

the world. First, the forest is felled and burned. Crops

are then planted in the ash and 2 to 4 harvests follow. The

land is then abandoned and crop production is shifted to a

new site. Shifting is necessitated by the marked, rapid

decline in crop production that almost invariably accompa-

nies cropping. The abandoned land is left in fallow for 4

to 20 years after which time the process is repeated.

The success of the system is based on the restoration

of soil productivity during the fallow period. If that

period is sufficiently long, the system works in an ecologic

sense (Sanchez, 1982). In much of the region, increasing

populations, coupled with governmental encouragement, is

causing influxes of peasant farmers and increased demand for

food. This, along with the desire to preserve some of this

unique habitat, is forcing reduced fallow periods and

consequent land degradation (Trenbath, 1984). These changes

are forcing the development and adoption of agricultural

techniques that are productive and can sustain continuous


Several factors have been attributed to the observed

yield declines with cropping. These include, declines in

soil fertility, increased weed infestation, deterioration of

soil physical properties, increased insect and disease

attacks, and social customs (Sanchez, 1976).

Whatever the cause of yield declines may be, it is

clear that, if sustained soil productivity is to be

achieved, soil fertility must be maintained or enhanced.

One way to achieve this is to add nutrients and other

required inputs in sufficient quantities to overcome yield

limitations. This approach has recently been adopted by

industrialized (mostly temperate zone) nations. Sanchez et

al. (1982) have shown that this approach can work in both an

agronomic and economic sense in a humid tropical setting.

The problem is that it requires a fairly well developed

infrastructure that provides farmers with access to a wide

variety of soil amendments, seed varieties and pesticides,

sufficient capital and credit to purchase them, sufficient

knowledge to properly manage them, and a well developed

marketing and transportation system.

An alternative approach that has received considerable

attention in recent years is based on the somewhat paradoxi-

cal observation that the land that supports some of the

world's most productive natural ecosystems (tropical

rainforest), supports some of the world's least productive

agro-ecosystems. If cropping systems were developed that

are more like the natural ecosystem, it is reasoned, they

should be more productive. This reasoning is supported by

the observation that most indigenous systems (excluding

paddy rice) that are of semi-permanent nature, are somewhat

like natural systems in that they are characterized by a

variety of different crops, including tree species, grown

together. The potential for designing such systems was

demonstrated by Hart (1980), who carefully planted crops to

mimic natural succession and obtained good, economic yields.

In practice it appears that few farmers are purists in

terms of the approach they adopt. A wide variety of crop-

ping systems and input levels are in use. In all cases, an

important consideration is that plant nutrients in the soil

be conserved. Low input systems will "run down" and high

input systems will be prohibitively expensive if nutrient

losses are high. One process that can cause nutrient losses

is leaching.

Leaching in Humid Tropical Regions

The observed soil productivity decline under shifting

agriculture has been attributed to leaching losses (Sanchez,

1976). Leaching is also considered to be a major factor

that limits the effectiveness of fertilizers in the region

(Engelstad and Russel, 1975).

At the outset it should be pointed out that the pro-

cesses involved in nutrient leaching in the tropical regions

are the same as those in temperate regions. A given nutri-

ent will be leached below the crop rooting zone when it

moves through the soil faster than the crop can absorb it.

This is dependent on the climate, soil and crop. Thus, the

amount of a nutrient entering the soil solution and the rate

at which it moves are important. Management practices, as

well as the soil, climate, and crop, effect this. Since all

of these factors are generally different in humid tropical

regions, the magnitude of the problem is potentially differ-

ent there than in temperate regions.

Rainfall is particularly noteworthy in this respect.

In many humid tropical regions the mean annual rainfall

exceeds 200 cm and, even though the potential

evapotranspiration is also high, it is often exceeded

considerably by rainfall. At Turrialba, Costa Rica, for

example, rainfall exceeds the potential evapotranspiration

by approximately 100 cm yr1 which is roughly equal to the

average annual precipitation in most temperate agricultural


Although generalizations about soils for such a large

region are dangerous, there are some differences between

soils of temperate and tropical regions that are common. In

particular, soils in humid tropical regions frequently

exhibit infiltration rates much greater than those expected

in temperate soils of similar clay content (Sanchez, 1976).

This is generally attributed to the relatively high degree

of aggregation commonly found in those soils. Since this

characteristic obviously affects water movement, it also

effects solute (nutrient) movement.

Another difference between soils is that the clay

fraction of the soils in the humid tropics is commonly

dominated by variable charged, low activity clays. This

effects the retention of both cations and anions.

Finally, the differences in individual crop species and

their management is frequently different from temperate

crops and management. All of these factors combine to make

transfer of agricultural experience and research in temper-

ate regions to the humid tropics tenuous. This is certainly

true of nutrient leaching studies.

Despite widespread concern that leaching may be an

important constraint to soil productivity in the humid

tropics, very little quantitative data have been reported in

the literature (Greenland, 1977; Keeney, 1982; Omoti et

al., 1983). This is partly because it is a difficult

parameter to measure and partly because much of the research

in the region is result-oriented (crop yield) as opposed to


Measurements that have been made range considerably.

Amounts of N estimated to have been lost in a year, for

example, range from 37 kg ha-1 (Omoti et al., 1983) under

oil palm, to 204 kg ha-1 from a pasture in Colombia

(Sanchez, 1976). A recent study of leaching from maize in

Nigeria calculated N losses to be about 80 kg ha- for a

single fertilizer application of 150 kg ha-1 and approxi-

mately half that for a split application (Arora et al.,

1982). Similar ranges of losses for K and Mg 2+, and even

greater ranges for Ca2+, have been reported (Sanchez, 1976).


At present it seems clear that there is a high poten-

tial for nutrient leaching losses in humid tropical regions.

It is also clear that research into management systems that

minimize such losses is an important priority. It has been

proposed that cropping systems can be manipulated in such a

way that such losses can be minimized. The principal objec-

tive of this study was to establish if nutrient losses from

two different cropping systems, which are strongly contrast-

ed in terms of species composition and diversity, are

different when levels of management are similar.

The fact that such wide ranges in leaching losses are

reported in the literature indicates that the processes

responsible for leaching must vary greatly in the region. A

second, related objective, was to quantify the processes

responsible for leaching in such a way that the sensitivity

of cropping systems can be characterized.


The two cropping systems chosen for comparison were:

(1) a mixed cropping system composed of laurel (Cordia

alliodora), cacao (Theobroma cacao) and platano (Musa

paradisiaca), and (2), a monocropping system composed of

maize (Zea mays L.). These two cropping systems will be

referred to as the MP, for mixed perennial, and MA, for

monocropped annual, in the remainder of the dissertation.

The two cropping systems were separated by approximate-

ly 100 m and were located in the "La Montana" section of the

experimental plots at CATIE (Centro Agronomico de

Investigaciones y Ensenanzas) near Turrialba, Costa Rica.

The soil series for both systems is Instituto clay loam

which is classified as a Typic Dystropept, fine, mixed,

isohyperthermic (Aguirre, 1971).

The management level of both systems was designed to

promote high yields and profitability. The MA plot is

located in a larger study area used by Dr. Carlos Burgos to

study the effects of tillage and residue management on maize

yields. It was initiated in November of 1976. Since then

40,000 plants ha-1 of a 120-day maize variety have received

approximately 400, 55, and 40 kg ha-1 of N, P, and K,

respectively, each year in four applications. Two crops

were grown each year. The dimensions of the plot discussed

in this study were 32 by 20 m.

The MP plot was part of a larger study of intercropping

with perennials directed by Dr. Gustavo Enriquez. It was

initiated in 1977. The planting densities were 1,111; 432;

and 123 plants per hectare for cacao, laurel, and platano,

respectively. The annual fertilizer regime was 140 kg ha-1

N, 30 kg ha-1 P, and 20 kg ha-1 K in four applications. The

dimensions of the MP plot were 18 by 18 m.


Three lines of inquiry, connected to the basic theme of

nutrient movement under the two cropping systems will be

presented in three separate chapters. First, investigations

of the nature of solute flow through the soil at the site

will be presented. This work was performed to establish the

appropriate model to be used in estimating field-scale water

and solute movement at the site.

The second line of inquiry is devoted to the develop-

ment and implementation of a soil-water and nutrient move-

ment model that was used to estimate leaching losses and

examine critical parameters associated with those losses.

The third line of inquiry is related to the supply of

nutrients to the systems via mineralization of soil organic

matter. In this section a laboratory method for estimation

of N mineralization is examined in terms of three models.

Implications of the results to field application of incuba-

tion results is also discussed.

The results presented in the three chapters are summa-

rized and discussed in terms of implications to the func-

tioning of the two cropping systems in the final chapter.



This study is part of a larger project investigating

nutrient movement under different cropping systems in Costa

Rica. Throughout the humid tropics there is a high poten-

tial for loss of plant nutrients by leaching due to the

considerable excess of precipitation over evaporative and

transpirational demand. Near Turrialba, Costa Rica, where

this study was performed, for example, the mean annual

precipitation of 264 cm exceeds mean annual pan evaporation

by 130 cm. Since improved crop yields are largely dependent

upon increasing the nutrient status of soils in the region,

quantification of leaching loss is important.

Quantitative description of solute (e.g., nutrient)

transport through packed soil columns has been accomplished

using the convective-dispersive (CD) model of solute flow

(Kirkham and Powers, 1972). Verification of this model has

come through numerous studies of miscible displacement

through uniformly packed, sieved, water saturated columns of

soil and other porous media (Nielsen and Biggar, 1961,

1962). According to the CD model, the transport of a

non-adsorped, conservative solute under steady water flow

conditions is described by

aC/at = Da2C/az2 v ac/az (2-1)

where C is the solute concentration (mg cm-3 ), z is distance

(cm), t is time (hr), D is the "effective" dispersion
2 -1
coefficient (cm hr ) and v0 is the average pore water

velocity (cm hr-1). The average pore-water velocity is

calculated by dividing the Darcy-flux (q) by the soil water

content 9. This implies that all soil water participates in

the convective transport of the solute.

In working with packed soil columns under steady water

flow conditions, two experimental conditions in which the CD

model has failed are (1) when aggregated media are used

(Nielsen and Biggar, 1961; Green et al., 1972), and (2) when

displacement was conducted under unsaturated conditions

(Biggar and Nielsen, 1962; Gupta et al., 1973; Gaudet et

al., 1977). In both cases this failure may be ascribed to a

failure of the assumption that all soil-water participates

in convective transport. When aggregated media are used,

intra-aggregate water is essentially stagnant so that virtu-

ally all convective transport occurs in the mobile water

located in the inter-aggregate portion of the soil. Simi-

larly, the application of tension to the soil causes the

entrance of gas which may result in the isolation of stag-

nant regions. Again, convective flux is restricted to the

mobile, nonstagnant-region water and Eq. 2-1 is not applica-


Application of miscible displacement techniques to

"undisturbed" soil columns taken from the field have shown

that inhomogeneities in field soils can strongly affect the

nature of water and solute flow (McMahon and Thomas, 1974;

Cassel et al., 1975). The presence of macropores is often

associated with solute flow behavior that is inconsistent

with the CD model (Kanchanasut et al., 1978). At present no

precise definition of macropores is widely accepted (Bouma,

1981; Luxmoore, 1981), but the term is generally used to

describe large, continuous pores that can conduct water and

solutes much more rapidly than the surrounding soil matrix

(Bevin and German, 1982; White, 1985). The process of flow

along macropores has been variously described as "channel-

ing", "bypassing", "short circuiting", "preferential flow"

or "partial displacement" (Scotter, 1978; Bouma, 1981; Beven

and German, 1982; Thomas and Phillips, 1979). The term

bypassing will be used in this paper.

The presence of relatively few macropores can greatly

increase the soil hydraulic conductivity (Bouma and Ander-

son, 1973) and cause very rapid transport of solutes through

the soil (Thomas and Phillips, 1979; Bouma et al., 1982)

that is inconsistent with the CD model (Kanchanasut et al.,

1978). The average pore water velocity (v ) is not an

appropriate descriptor of convective transport when there is

bypassing because flux is effectively limited to a small

portion of the total soil-water (i.e., that within


Application of Eq. 2-1 on a field-plot scale has shown

that the parameters v0 and D are extremely variable and

log-normally distributed (Biggar and Nielsen, 1976), which

makes their estimation very difficult. Application of

numerical models, which have been developed for more realis-

tic, transient conditions that exist in the field, are

subject to similar limitations.

As an alternative, management level models that use

less variable, capacity-type parameters (e.g., 8) have

proven successful in many cases (Rao et al., 1981; Rose et

al., 1982). These models are largely empirical in nature

but utilize the concept of convective transport, as defined

by v Failure of these simplified management level models

has been attributed to flow along macropores (Thomas and

Phillips, 1979; Barry et al., 1985).

The three situations mentioned above in which the

inapplicability of Eq. 2-1 have been documented; unsaturated

flow; flow around aggregates; and flow along macropores,

have in common the fact that a portion of the soil-water,

called mobile water (8m ) in this paper, moves much more

rapidly through the soil than the remaining, immobile water

(9 im). Surface applied solutes can be transported much more

rapidly in the mobile region and thereby bypass solutes

residing in immobile regions.

An alternative model, known as the mobile-immobile

water model (MIM) has been developed to explicitly account

for this situation (Passioura, 1971; Van Genuchten and

Wierenga, 1976). In this model two soil-water phases are

explicitly differentiated, with all convective transport

assumed to occur in the mobile phase and transport in the

immobile phase restricted to diffusion. Since no specific

pore geometry is assumed, the diffusive transfer between the

two soil-water phases is described as being proportional to

the concentration difference between them. The equation for

transport of a nonadsorped solute under steady water flow

conditions that is consistent with the MIM model is

8 mC/at + 9 imCi /At= D2C /az2 8 mv9C/az (2-2)

with interphase transfer described as

9im Cim/at = a(C -C im) (2-3)
where the subscripts "m" and "im" refer to mobile and

immobile phases, respectively. The parameter a is an

empirical constant called the mass transfer coefficient
(hr~ ). The parameter vm differs from vo in that it is

calculated from q/8 Simplified, capacity-parameter based

management level models, analogous to those based on the CD

model, have been developed for the MIM model (Addiscott,

1977). Use of this model requires the fraction of mobile

water 0 (em/9), as an input.

Our primary objective in this study was to select a

modeling approach to describe water and solute movement

under field conditions. From the above discussion it is

clear that the nature of the conducting pore network, or the

hydrologically effective pore geometry, greatly affects the

applicability of a given model. Ideally, it would be

possible to infer which model (if either) is appropriate

from standard soil descriptions. Although the impact of

soil structure on the hydrologically effective pore geometry

is well documented (Elrick and French, 1966; Cassel et al.,

1974; McMahon and Thomas, 1974), and the effects of soil

texture are widely recognized (Thomas and Phillips, 1979;

Addiscott and Wagenet, 1985), the information is of a

qualitative nature (Bouma, 1981). More detailed morphologi-

cal descriptions are generally not well suited to this

purpose because they do not yield information on pore

continuity or the degree of pore interconnections (Bouma,

1981). For these reasons, study of the hydrologically

effective pore geometry has usually focused on measurements

of soil hydraulic properties (Bouma and Anderson, 1973) and

movement of tracers (Bouma, 1981; White, 1985).

Our secondary objective was to identify soil structural

or pore characteristics that characterize the hydrologically

effective pore geometry. This information is potentially

useful in using field-level observations to infer the nature

of flow through soils. Dyes have been used extensively for

this purpose (Bouma and Dekker, 1978; Omoti and Wild, 1979).

Materials and Methods

Soil Characteristics and Management

The soil studied is an Instituto clay loam (fine-loamy,

mixed, isohyperthermic, Typic Dystropept) located at CATIE

(Centro Agricola Tropical de Investigaciones y Ensenanzas)

near Turrialba, Costa Rica. It is derived from alluvium

deposited from the surrounding mountains which are primarily

of volcanic origin. Very little soil profile development

was evident except for an accumulation of organic matter at

the surface. Textures are clay loam throughout. Aguirre

(1971) described the structure as weak, subangular blocky

with peds ranging from 0.5 to 2.0 cm in diameter near the

surface and becoming increasingly less pronounced with

depth. The site is nearly level and the soil is considered

to be moderately well drained (Aguirre, 1971).

Soil columns were taken from two experimental plots

located about 100 m apart. One plot had been planted to

maize for six years prior to column removal. The crops on

the other plot were a mixture of cacao, laurel, and plantain

and had been under continuous management for five years

prior to column removal. No tillage was performed on either

plot during that time and no machinery entered either plot.

Columns numbered 1, 2, and 3 (Table 2-1) were taken from the

first plot and those numbered 4 and 5 were taken from the


Column Experiments

Five undisturbed soil columns were collected at two

depths (0-30 cm and 75-105 cm). The procedure used to

collect the samples was as follows: (1) dig a pit approxi-

mately 1.0 by 1.5 m leaving a pedestal approximately 0.3 m2

in the middle, (2) cut a circle about 15 cm in diameter in

the middle of the pedestal with a sharp knife, (3) force the

12 cm diameter PVC tube that is beveled at the end and lined

with petroleum jelly into the area bounded by the circle,

(4) shave the surrounding soil from the edges of the pipe

and cut another circle about 1 cm deep, (5) repeat the last

two steps until the 30 cm long column is filled, and (6)

separate the column from the soil at the bottom. The upper

12 to 15 cm of each column were used in the displacement


Fritted glass endplates were fitted to both the inlet

and outlet ends of the columns. A constant hydraulic

potential was maintained at both column ends during dis-

placement experiments. Tension was applied at the inlet end

via a Marriot devise and at the outlet end by a hanging

water column. All displacement experiments were conducted

under a unit hydraulic head gradient except the two dis-

placements performed in Column 1 (Table 2-1) under saturated

conditions. The influent 3H activity was about 7 nCi ml-1

in a 0.01 M CaCl2 solution. The 3H activity in effluent

fractions was assayed using liquid scintillation techniques.

Rhodamine B dye displacements were performed using the

same columns after the 3H20 displacements. The influent dye

concentration was 2 g L 1. The amount of dye displaced

ranged from 0.05 to 0.32 pore volumes (Table 2-1). Dye

patterns in cross-sections of the columns were photographed

at 1 cm intervals after displacement. Selected

cross-sections were then traced with pen and ink.

The experimental conditions under which all displace-

ments were performed are summarized in Table 2-1.

Table 2-1. Experimental conditions for displacement.

Column Column Depth Tension of Tension of Pore Vol.
Number Length H Displ. Dye Displ. Displ.*
cm ---------kPa-----------
1 12.5 surface 0.0,0.1,0.5, 1.0 0.32
2 15.0 subsoil 0.0,1.0 0.0 0.05
3 12.5 surface 0.0,1.0 0.0 0.13
4 12.0 subsoil 0.0,1.0 1.0 0.19
5 12.5 surface 0.0,1.0 ---
* The number of pore volumes displaced during the dye

Field Techniques

The depth to the water table was measured in perfo-

rated, 2.54 cm diameter PVC pipe. Measurements were made

several times each week. Rainfall was measured daily at the

site. Infiltration measurements, carried out under the

direction of Dr. Carlos Burgos, were made with a double-ring

infiltrometer at 24 locations on an adjacent plot. The

dimensions of the inner and outer rings were 31 and 60 cm,


Adsorption Experiments

Adsorption of 3H20 was measured using two different

batch techniques under two slightly different conditions.

The first method was that described by Dao and Lavy (1978)

in which the soil-solution ratio was 2 to 1 (g g- ). The

soil used was taken from column 3 after the displacement

experiments had been performed and the soil was oven dried

for pore volume determination.

In the second method, air-dried soil from column 3 but

not used in the displacement experiment, was added to 3H20

solution in a soil-solution ratio of 1 to 2 and placed in a

shaker overnight. Concentrations of 3H20 ranged from 7 to

0.7 nCi ml-1 in both experiments.

Parameter Estimation and Model Evaluation

In Eqs. 2-1 through 2-3 no adsorption was assumed,

which is unrealistic for most tracers. Expansion of Eq. 2-1

to include adsorption results in the following expression:

DC/;t + (BD/8)(OS/at) = D92C/Iz2 v aC/az (2-4)

where BD is the bulk density (g/cm 3), S is the adsorbed

phase concentration (mg g- ) and the other parameters have

been previously defined.

Adsorption can frequently be described by a linear (or

linearized) equation of the form

S = KDCe (2-5)

where Ce is the equilibrium solution concentration (mg cm-3)

and KD is an empirical distribution coefficient (g cm-3 ).

If the adsorption process is assumed to be instantaneous and

reversible, Eq. 2-5 can be substituted into Eq. 2-4 to give

RFaC/At =DP2C/z2 -v 0C/8z (2-6)

where RF, the retardation factor, is defined by

RF = 1 + BDKD/8. (2-7)

Analysis and interpretation of parameter values is

facilitated by the use of dimensionless variables. Equation

2-6 can be described in terms of the following dimensionless


T = v t/L, (2-8)

x = z/L, (2-9)

P = v L/D, and (2-10)

C = Cb/C, (2-11)

where v t, z and D have been defined previously, L is the

column length, P is the Peclet number, T is the number of

pore volumes, x the dimensionless distance, and C the ratio

of effluent (Cb) to influent (C ) concentration. (Note that

the initial soil solution concentration is assumed to be 0).

Substitution of these variables into Eq. 2-6 results in the

following expression:

RF(aC/3T) = (l/P)(a2C/aX2) aC/ax. (2-12)

In order to incorporate adsorption into Eqs. 2-2 and

2-3, S is partitioned between mobile and immobile phases

such that

S = fSm + (1-f)Sim (2-13)

where Sm and Sim are the adsorbed concentrations in the

mobile and immobile regions, respectively, and f is the

fraction of adsorption sites in the mobile region. Assuming

that the same linear relationship expressed in Eq. 2-5

applies to both the mobile and immobile regions, Eqs. 2-2

and 2-3 can be written

(9m+BDfKD)aCm/9t + [im + (l-f)BDKD]'Cim/3t =

(86D)32Cm/3x2 (9mv )3Cm/ax (2-14)

[9i+(l-f)RoKD]aCim/3t = a(Cm-Cim). (2-15)
As with the CD model, the MIM model can be described in

terms of dimensionless variables. The variables used are

P = (9m + fBDKD)/(9 + BDKD), (2-16)
w = aL/q, (2-17)

c1 = Cm/C', and (2-18)

C2 = Cim/Co, (2-19)
where cm and Cim refer to soil solution concentrations in

the mobile and immobile regions. All other parameters have

been defined previously. The MIM model for steady-state

water flow conditions can now be described by

3RF(Dcl/3T) + (1-3)RFacl/DT =(1/P)(D2cl/ax2) acl/ax (2-20)


(1-3)RFac2//T = w(c1-c2). (2-21)

The result of these transformations is that the CD

model is described by Eq. 2-12 and the MIM model by Eqs.

2-20 and 2-21. Measured solute breakthrough curves (BTC's)

were fit to both the CD and MIM models. The program CFITIM

(van Genuchten, 1981), which uses a nonlinear, least sum of

squares criteria for goodness-of-fit was used. In every

case a first type, constant concentration, influent end

boundary condition was assumed. The other boundary condi-

tion was that of a semi-infinite column.

Results and Discussion

Hydraulic Parameters

The Darcy flux (q, cm hr-1), soil-water content (89, cm3
cm 3), and soil-water tension (h, kPa) under which the

experiments were performed are shown in Figs. 2-1 through

2-5. The hydraulic conductivity (K, cm hr 1) value is

included for the two BTC's in which K was not equal to q.


Hydraulic conductivity values were extremely high

considering that the soil has a clay content of greater than

30%. The combination of high clay content and high saturat-

ed hydraulic conductivity (Ksat) is frequently indicative of

flow along macropores (Bouma and Anderson, 1973; McKeague et

al., 1982). Saturated 9 values are also fairly high.

A dramatic decrease (two orders of magnitude) in K was

observed as h increased from 0 to 2 kPa (Fig. 2-1). This

was accompanied by a relatively modest decrease (6.5%) in 8.

Similar trends in K and 8 were observed in the other columns

(Figs. 2-2 to 2-5). This indicates either that: (1) few

discrete, large (>0.5 mm radius) pores conduct large volumes

of water rapidly under saturated conditions, or (2) water

held in large pores serves to "connect" a number of pores

that conduct water rapidly.

The series of tensions in Fig. 2-1 illustrates the

transition between saturated and unsaturated conditions for

the surface soil. Preliminary investigations with subsoil

showed a qualitatively similar pattern. This similarity

between subsoil and topsoil was confirmed in subsequent

comparisons between saturated and unsaturated conditions

(Figs. 2-2 through 2-5).

Qualitative Evaluation of Breakthrough Curves

In general, BTC shapes changed dramatically as soil

water tension was increased. BTC's obtained under saturated

conditions were highly skewed, characterized by very early


.a .o o b ..--)" C

C/C h= 0.0 h = 0.0
o 0= 0.57 8=0.57
q = 26.4 q = 2.7
0.2 K= 40.2 K= 40.2

c d

.C/Co h= 0. I h= 0.5
09=0.57 9=0.55
q = 1=9.5 q = 5.9

1.0 e 0

C/C h= 1.0 h= 2.0
0= 0.53 0= 0.52
Sq= 1.2 = 0.2

1.0 2.0 3.0 1.0 2.0 3.0

Figure 2-1, a-f. The effect of soil-water tension (h, kPa),
soi}-water content (9, cm cm ), and Darcy flux (q, cm
hr ) on the elution of tritiated water in Column 1. The
continuous line represents the best fit of the MIM model in
a, b, and c, and of the CD model in d, e, and f.

1.0 b-- oO

C/C h=O .0 h= =l.0
S -=0.60 8= 0.59
Sq= 28.7 Iq= 1.3
I f
1.0 2.0 3.0 1.0 2.0 3.0

Figure 2-2, a and b. Effect of soil-water tension (h, kpa),
soil-water content (9, cm cm ), and Darcy flux (cm hr )
on elution of tritiated water in Column 2. The solid lines
in a and b represent the best fit of the MIM and CD models,


1.0 2.0 3.0 1.0 2.0

Figure 2-3, a and b. The effect5 of soil-water tension (h,
kPa),_ oil-water content (9, cm cm ), and Darcy flux (q,
cm hr ) on elution of tritiated water in Column 3. The
solid lines in a and b represent the best fit of the MIM and
CD models, respectively.

o 9~8=0.54 9=0.52
q= 29.4 q= 0.9
0.2 ?

1.0 2.0 3.0 1.0 2.0 3.0

Figure 2-4, a and b. The effect of soil-water tension (h,
kPa),_ oil-water content (9, cm cm ), and Darcy flux (q,
cm hr ) on elution of tritiated water in Column 4. The
solid lines represent the best fit of the MIM and CD models
in a and b, respectively.


.0 2.0 3.0 1.0 2.0

Figure 2-5, a and b. The effect of soil-water tension (h,
kPa),_ oil-water content (9, cm cm ), and Darcy flux (q,
cm hr ) on the elution of tritiated water in Column 5. The
solid lines represent the best fit of the MIM and CD models
in a and b, respectively.

appearance of tracer in the effluent and a slow approach of

effluent concentration towards 1.0 (also called "tailing").

This behavior has been observed in undisturbed columns

(Anderson and Bouma, 1977; White et al., 1984), and has been

inferred from field studies of solute movement (Wild and

Babiker, 1976). Such skewed BTC's indicate that solute was

conducted relatively rapidly through the columns by a small

fraction of the total soil-water. The rapidly conducting

fraction of the soil-water has been related to

inter-aggregate regions (Nkedi-Kizza et al., 1982),

inter-ped regions (Anderson and Bouma, 1977), and discrete

macropores (Kanchanasut et al., 1978) in other soils.

Unsaturated BTC's are markedly more symmetric, with a

later arrival of tracer and less tailing. This trend was

not reversed as tension was increased to 2 kPa. These

results contradict other studies that have shown that

tailing in BTC's resulted from increases in soil water

tension (Nielsen and Biggar, 1961; Gaudet et al., 1977).

Apparently the pore network in this soil is sufficiently

interconnected that drainage of large pores does not result

in the isolation of stagnant regions in the column. It

should be noted that both of the studies referenced above

were performed on sands in packed columns at relatively low

soil-water contents. Lower soil-water contents enhance the

possibility of the creation of isolated, stagnant regions in

the soil. One study by Elrick and French (1966) that

compared saturated and unsaturated flow in an undisturbed

column found that dispersion decreased with application of

tension, although marked asymmetry during saturated flow was

not observed.

The observed change in BTC shape with increasing soil

water tension can be explained by the concomitant decrease

in q and/or by changes in the effective pore geometry. In

general, bypassing is enhanced by greater fluxes at a given

8 because because there is less time for diffusive transfer

into stagnant regions. At the same time, increasing tension

changes the effective pore geometry by draining larger pores

that may be responsible for the observed bypassing.

Both effects were operative in this study but the

effect on the effective pore geometry was dominant. This is

illustrated in Fig. 2-1. When flow rates were reduced and

the soil remained saturated (Figs. 2-la and b), the BTC

shape was only slightly altered. However, when displacement

on the same column was performed under unsaturated condi-

tions at approximately the same q (Figs. 2-lb versus 2-ld

and e), there was a considerable change in BTC shape. These

trends were reflected in the models and parameters used to

describe the BTC's.

Quantitative Evaluation of Breakthrough Curves

All BTC's were fit to both the MIM model and the CD

models. In general, curve fits fell into two groups, satu-

rated and unsaturated, with the 0.1 kPa run (Fig. 2-lc)

intermediate. The parameter values obtained will be dis-

cussed by these groups.

Unsaturated Breakthrough Curves

Two dimensionless parameters are required in the CD

model, P and RF. The best least sum of squares fit was

obtained allowing both parameters to vary. The solid line

in all the unsaturated BTC's except the 0.1 kPa run repre-

sent the calculated best fit using the CD model (Figs. 2-1

through 2-5). In general, the agreement was excellent. The

resultant parameter values and associated 95% confidence

intervals shown in Table 2-2 appear to be independent of the

depth from which the columns were sampled. Peclet numbers

(P) ranged from 3 to 12 and RF values from 1.12 to 1.17.

Dispersion coefficients calculated from those P values were

high relative to those measured in sieved, packed columns,

but this is expected in undisturbed columns (McMahon and

Thomas, 1974; Cassel et al., 1975). It indicates that there

was a relatively wide range of pore water velocities within

the column.

The RF values obtained were high considering that

tritium is frequently assumed to be nonadsorped (RF=1.0).

Tritium sorption has been noted by several workers (Mansell

et al., 1973; Wierenga et al., 1975; Van de Pol et al.,

1977; Nkedi-Kizza et al., 1982) and has been associated with

hydroxyl exchange with clay lattice hydroxyls (Stewart,


As a check of the accuracy of RF values determined by

curve-fitting, the area above the measured BTC was calculat-

ed for several of the unsaturated BTC's (Pandey and Gupta,

1984). The areas measured agreed closely with the RF values

obtained by curve-fitting. Since the fit between calculated

and measured curves was excellent in every case, this result

confirms the finding of van Genuchten and Parker (1984) that

mass balance is preserved with the solution and boundary

conditions used.

Table 2-2. Convective-dispersive model parameter values.

Col h P RF KD D
-1 -
kPa ml g- cm2 hr1
1 0.1 0.90 1.08 0.031 492
(0.02)* (0.01) (0.006) (11.7)
1 0.5 2.28 1.10 0.043 60.9
(0.06) (0.01) (0.006) (0.02)
1 1.0 4.45 1.12 0.053 6.44
(0.31) (0.02) (0.008) (0.17)
1 2.0 6.62 1.17 0.069 0.58
(0.17) (0.01) (0.002) (0.01)
2 1.0 13.5 1.12 0.062 2.52
(0.76) (0.01) (0.005) (0.13)
3 1.0 13.1 1.18 0.080 3.78
(0.49) (0.01) (0.002) (0.14)
4 1.0 12.2 1.16 0.071 1.65
(0.57) (0.01) (0.003) (0.07)
5 1.0 7.01 1.18 0.079 5.81
(0.36) (0.01) (0.005) (0.18)
* Numbers in parenthesis are the 95% confidence intervals
associated with the estimated value.

An independent check of RF can be obtained from mea-

surement of adsorption in batch isotherms. The batch

isotherms obtained using both methods described in the

previous section were essentially identical. Both isotherms

were linear with r2 values of 0.994 and 0.997 and KD values,

with 95% confidence intervals of 0.134 + 0.0045 and 0.132 +

0.0047 mg g-1. These values are significantly larger than

those obtained from values derived from fitted parameters

(Table 2-2).

Discrepancies between batch and column-measured adsorp-

tion parameters have been noted by others (Nkedi-Kizza et

al., 1982). The value of RF calculated using the

batch-derived KD value (including a coarse fragment content

of 3.9% in the column) is 1.26, which is slightly higher

than that obtained from curve fitting (1.17). This discrep-

ancy is likely due to differences in the condition of the

soil when the experiments were performed. The batch iso-

therms were conducted on air-dried or oven-dried soil while

the columns were never air-dried.

When all unsaturated BTC's except the 0.1 kPa run were

fit to the MIM model either extremely high w values (>35) or

3 values of 1 were obtained. From inspection, it is clear

that Eqs. 2-20 and 2-21 are indistinguishable from Eq. 2-12

when P=l, and use of the MIM model is not justified. The

effect of high w values is less clear but implies extremely

rapid transfer between mobile and immobile regions which

effectively eliminates the need to make a distinction

between them. This will be discussed in greater detail in

the next section.

Saturated Breakthrough Curves

The 0.1 kPa run (Fig. 2-1c) will be included in this

part of the discussion because it more closely resembles the

saturated BTC's than the other unsaturated BTC's. When the

CD model was applied to the saturated BTC's generally poor

fits resulted. This was due to the very rapid rise in C/C0

and obvious leftward shifting. Use of the MIM model

requires specification of 4 dimensionless variables; P, RF,

3, and w. Although it is possible to allow all variables to

vary simultaneously, the resultant parameter estimation is

relatively imprecise. One value, RF, should be consistent

with those derived from the unsaturated (CD model) BTC's.

This value was accordingly taken from the unsaturated curves

and fixed during parameter estimation for the saturated


In every case, very close agreement between measured

and calculated BTC's was obtained. In general, the parame-

ter values in Table 2-2 indicate the following trends; very

small P values corresponding to extremely large D values; 3

values on the order of 0.23 to 0.45; and w values ranging

from 0.33 to 3.84.

Table 2-3. Mobile-immobile water model parameter values.

Col. P RF P w D ESR D

cm2/hr cm
1 (fast) 0.29 1.12 0.76 1.02 2515 0.18 0.81
(0.01)* (0.03) (0.39)
1 (slow) 0.86 1.12 0.43 2.46 150 0.61 0.47
(0.05) (0.04) (0.45)
1 (1kPa) 1.11 1.12 0.84 0.08 437 0.51 0.89
(0.03) (0.04) (0.03)
2 0.43 1.13 0.34 3.84 1643 0.18 0.37
(0.03) (0.03) (0.48)
3 0.55 1.17 0.34 2.99 2758 0.19 0.38
(0.03) (0.03) (0.48)
4 0.76 1.15 0.23 0.16 3497 0.86 0.26
(0.30) (0.05) (0.03)
5 0.14 1.17 0.45 3.89 7525 0.15 0.51
(0.02) (0.11) (2.53)
* Numbers in parenthesis are the 95% confidence intervals
associated with the estimated value.

The very low P values are indicative of an extremely

broad range in pore-water velocities in the mobile water

region. This indicates that the compartmentalization of

soil-water into only two phases assumed in the MIM model was

insufficient to account for the range of pore-water veloci-

ties encountered. It is possible that the logic of the MIM

model be extended to consider gradations of soil-water

mobility (e.g., "very rapidly mobile water", "somewhat

rapidly mobile water", etc.). Skopp et al. (1981) have

applied this approach but to only two soil-water phases.

Recent work applying the transfer function model developed

for soil applications by Jury (1982) to the analysis of

BTC's similar to the saturated BTC's presented here (White

et al., 1986) have described the distribution of pore-water

velocity as a continuous, albeit skewed, function. Viewed

in this perspective the MIM model may be considered to

represent an extreme, bimodal pore-water velocity distribu-


Calculation of D from fitted P values is somewhat

questionable when P values are so low. There are two sets

of boundary conditions that approximate the experimental

conditions and the solutions to those boundary conditions

diverge when P is less than 4 or 5 (see van Genuchten and

Parker, 1984; Parlange et al., 1985).

The parameter P means little on its own, but can be

related to the mobile water fraction, 0 (recall that 0=9 /8)

with the expression 0=(RFB)-f(RFl). Note that, if RF=1,

3=P3. As a first approximation, 3 can be considered to be a

measure of $ in these experiments because RF is close to 1.

A more refined estimate of Q is obtained if some assumptions

concerning f are made. Recall that f was defined as the

fraction of sorption sites in the mobile region. If the

distribution of sorption sites is independent of location in

soil-water regions, then f=O when the soil is saturated and

D=A. However, it seems reasonable to expect that propor-

tionately more sorption sites will be found in immobile than

mobile regions because the pores in immobile regions should

be smaller and therefore have more exposed surface area.

This reasoning has been used to justify the assumption that

f=O (Nkedi-Kizza et al., 1982) which leads to O=RF3. Thus,

in soils with positive adsorption, $ values of greater than

0 are expected. The 0 values in Table 2-3 were calculated

assuming that f=0/2, which is an intermediate estimate.

The results in Table 2-3 indicate that $ values of 0.25

to 0.50 are generally consistent with the parameters fitted

with the MIM model. These values are surprisingly high in

light of the large changes in K that resulted from relative-

ly small changes in 9. An independent estimate of 0 can be

obtained by assuming that water held at field capacity is

immobile (Addiscott et al., 1977). In this case, such an

approach yields estimates of 0 of approximately 0.156, which

is distinctly lower than those obtained from curve-fitting.

When this value was fixed along with RF the resultant P

values were increased, w values were decreased and the

goodness-of-fit was substantially reduced. These fitted

curves were too angular, displaying a more rapid rise in

effluent concentration with greater tailing than the mea-

sured BTC's. The difficulty of obtaining independent

estimates of $ has been noted by others (Addiscott et al.,


The parameter w is somewhat more difficult to interpret

as it is not directly related to any specific soil charac-

teristic or property. However, work by Rao et al. (1980a)

has shown that w can be used to successfully calculate

inter-aggregate concentrations during diffusion into spheri-

cal aggregates of known volume. In subsequent work (Rao et

al., 1982) it was demonstrated that media composed of

different sizes and shapes of aggregates could be approxi-

mated by a single "equivalent" spherical aggregate size on a

volume-weighted basis. This work has recently been extended

to miscible displacement studies by van Genuchten (1985).

Using analytical solutions for flow through media

composed of immobile regions of known geometry van Genuchten

(1985) was able to express w in terms of an average sphere

or other aggregate shape. This technique was applied to the

saturated BTC's to determine size of effective spherical

radius (ESR) consistent with the fitted parameter values.

The parameters shown in Table 2-3 indicate that the soil may

be considered to be composed of spherical aggregates of of

0.15 to 0.6 cm radius.

Another way of considering the effect of aggregate size

on the effective pore geometry was presented by Rao et al.

(1980b). They showed that, when the aggregate size is small

enough relative to the pore-water velocity, a condition of

near-equilibrium will be established and the CD model should

be appropriate. For spherical aggregates, the condition of

near-equilibrium is valid when

DeL(1-)/(a2v 020.3) > 1 (2-22)

where De is the diffusion coefficient and L is the column

length (Rao et al., 1980b). Taking ESR as 0.5 leads to a

critical v0 of 1.6 cm/hr, which is generally less than the

v values of the unsaturated runs (v =2q in these columns).

Thus spherical regions of immobile water could exist in the

columns but their effect on solute transport would be

"masked" as high dispersion.

Taken as a whole, some inferences concerning the nature

of the effective pore geometry of the soil can be made. The

high Ksat values measured in conjunction with highly skewed

BTC's indicate that water was conducted relatively rapidly

through some portion of the soil when saturated. The large

reduction in K and BTC skewedness that resulted from appli-

cation of 0.1 to 0.5 kPa of tension suggests that soil water

was rapidly conducted via macropores (Luxmoore, 1981;

Germann and Bevin, 1981). However, the extremely low P

values and 0 values of 0.25 to 0.50 under saturated condi-

tions suggest that there were several such regions of

varying conductivities. The model parameter values obtained

are consistent with fairly large regions of immobile water.

If immobile regions are assumed, for example, the immobile

regions would have radii of 0.15 to 0.6 cm.

Although immobile regions were described in terms of

spherical aggregates, the MIM model specifies no pore

geometry and numerous other possibilities exist. Rhodamine

B dye was used to better determine the actual effective pore


Dye Experiment

Dyes have frequently been used to visually investigate

the nature of flow paths through the soil (Bouma and Dekker,

1978; Omoti and Wild, 1979; McVoy, 1985). The approach

provides the opportunity to observe the conducting pathways

and thereby relate water and solute flow to observable

structural features or biochannels. The basic assumption

made in interpreting dye patterns is that the more solution

that passes a given point, the more darkly stained that

point will be. Thus, stained regions are interpreted as

being regions of relatively fast flow, unstained regions to

be of relatively slow flow.

It is important that this fairly simple-minded approach

not be extrapolated far in terms of correlation of dye

patterns with BTC's. In the first place, Rhodamine B dye is

sorped to the soil much more strongly than tritium (McVoy,

1985) so that dispersion is apparently reduced. Secondly,

the dye is not instantaneously and reversibly desorped as

tritium is assumed to be. And thirdly, visual evaluation of

color is qualitative, so that quantification of the amount

of dye at a location is not possible.

Given these difficulties, four observations of note can

be made from the dye patterns illustrated in Figs. 2-6

through 2-9. First, no significant staining of the column

edges was observed in Columns 1, 2, and 3 and the staining

on the edge of Column 4 was not as intense as in the inter-

nal portions of the column. This is a critical question

that must be considered when undisturbed columns are used.

From these observations we do not believe that observed

BTC's were strongly affected by the presence of the column


Second, while specific stained regions that must have

been responsible for the very early appearance of tracer in

the effluent were easily identified, with one exception,

they were not obviously associated with discrete biochannels

or structural features. Even with the segmented column at

hand, it was very difficult to determine exactly which pores

were conducting, as, in every case many visible pores were

stained. In addition, it was difficult to trace individual

pore sequences up the column because they meandered consid-

erably across the column. Omoti and Wild (1979) and McVoy

(1985) have made a similar observations.

Third, where structural units were relatively strong as

in Column 1, there was preferential flow around them. The

structural units isolated in Column 1 ranged in radius from

about 0.3 to 1.2 cm, which is in rough agreement with that

0 --

n,A W

9- -
12 -


0 5 10 15

Figure 2-6. Rhodamine B dye staining pattern in Column 1
resulting from displacement under a soil-water tension of 1


=s -- -

/- ----~

% -O.'
=r --

-. -. .. .0. _- .-

-- I J
-. I,C

4^ f ^

3 5 10 15

Figure 2-7. Rhodamine B dye staining pattern in Column 2
resulting from displacement under saturated conditions.




Figure 2-8. Rhodamine B dye staining pattern in Column 3
resulting from displacement under saturated conditions.

*i -.m y *%

- .. '

r A ys
i ^" ?^ '
' ^e~k .
-- *-

-~ -

- ---

o 5 10 1!

Figure 2-9. Rhodamine B dye staining pattern in Column 4
resulting from displacement under a soil-water tension of 1






calculated from the curve-fit parameters. The fact that

these units were identified during unsaturated flow which

was well described by the CD model indicates that movement

into and out of those units was sufficiently rapid that

bypassing was not indicated in the BTC. This observation

probably accounts for the generally greater dispersion

observed in undisturbed columns and field studies. That is,

bypassing is "masked" in the dispersion term.

Fourth, the nature of the stained regions does not

appear much different in the saturated and unsaturated

columns. This is evidence that, rather than drain a few,

discrete large pores, the application of tension drains

regions of the soil that serve to "connect" pore-sequences.

The main difference between saturated and unsaturated

dye patterns is that the solutee front" is more compressed

in the unsaturated columns. Note that the number of pore

volumes of dye solution applied to all columns was approxi-

mately the same (Table 2-1), but the extent of staining in

the unsaturated runs was more strongly weighted toward the

inlet end of the column.

These observations compliment the results of BTC

analysis and hydraulic property measurement in the previous

section. It appears that the highly skewed BTC's and very

high Ksat measured under saturated conditions were due to

very rapid transport in a very restricted region of the

columns. These regions are better characterized as conduct-

ing pore sequences than discrete macropores and their

identification in the field would be very difficult. There

appears to be a number of such pore sequences that range

widely in conductivity. Application of tension "discon-

nects" the largest effective pore sequences and therefore

results in reduced skewing and K. Immobile regions were

generally characterized as regions between conducting pore

sequences as opposed to easily identifiable, physically

controlled regions. When displacement occurred under

tension, flow in the conducting regions was slow enough, and

those regions were close enough, that observed heterogeneous

flow was described as dispersion with the CD model.

Field Application

One of the stated objectives of this study was to use

information from column experiments as a basis for selecting

a model to describe movement of nutrients and water under

field conditions. The basic distinction between the two

models considered is whether or not all soil-water effec-

tively participates in convective transport (i.e., whether

or not v0 applies). It is clear from the dye patterns in

Figs. 2-6 through 2-9 that there was bypassing in the sense

of heterogeneous flow under both saturated and unsaturated

conditions. In terms of model selection, however, bypassing

was significant only when the soil-water content was at or

very near saturation. Hence, "significant" bypassing is

expected in the field only when the soil is near saturation.

The necessary precondition for saturation of the soil

surface under field conditions is that the water input rate

(rainfall intensity) exceed the infiltrability of the soil

(Hillel, 1971). If no surface disturbance occurs (none was

observed) the minimum infiltrability of a uniform soil

profile is K sat. Since no differentiation of horizons in

terms of hydraulic properties was observed, the rainfall

intensity must at least exceed K sat if saturation and hence

"significant" bypassing are to occur.

The Ksat values measured in undisturbed columns were

between 20 and 40 cm hr-1. Given the well known high

spatial variability of Ksat (Warrick and Nielsen, 1980),

these values probably should not be used in this context.

We obtained estimates of Ksat from 24 double ring infiltra-
tion measurements. The results are shown in Fig. 2-10. The

values obtained are considerably lower than those measured

in the columns. Aside from spatial variability, two

explanations for this difference are the fact that the

columns were saturated from below under tension and were

therefore closer to true saturation and the fact that

macropore continuity is enhanced in shortened columns

(Edwards et al., 1979)

When the frequency distributions of Ksat and rainfall

intensity are compared (Fig. 2-10), it is evident that

rainfall intensity sufficient to cause ponding is rare.

Based on this information it appears that saturated flow and

the attendant bypassing are not expected to occur at the

site. This is consistent with the observation at the site

that ponding did not occur.

I I\

_D ^




4 -



(suoieBAJasqo %) AuenbaJj OA4ieJOa


aC c

4J 0)

-I .



*H 4


c a
r) -H



N *U
.H (
3 f








Application of BTC's measured in columns to field

conditions involves a considerable extrapolation in scale.

This requires that the representative elementary volume

(Bear, 1972) for bypassing be considerably less that the

volume of the column. Several studies have shown that

measured bypassing at a field scale was expressed in undis-

turbed soil columns (Bouma and Wosten, 1979; Omoti and Wild,

1979; White 1985). This, however, need not be the case.

If significant bypassing at a scale larger than column

dimension is to occur, there must be some means by which

soil water movement is concentrated. This may be the result

of soil properties or the distribution of incoming water.

Large soil structural units or biochannels have been shown

to be responsible for bypassing (Bouma and Dekker, 1978;

Bouma et al., 1982). These effects are enhanced by a very

slow matrix Ksat. Other soil characteristics that could

serve to concentrate flow are coarse fragments, steep

slopes, and strongly contrasting horizons (Bevin and

Germann, 1982).

The soils in this study exhibited none of the above

characteristics. Soil structural units were generally less

than 2 cm in diameter and poorly expressed. From the column

studies it is clear that the matrix conductivity is rela-

tively high (>1 cm hr-1). Other characteristics such as

steep slopes, contrasting horizons and coarse fragments were

not evident. Soil animals that could potentially make large

channels (e.g., leaf cutter ants and armadillos) were

carefully excluded from the site.

The input of water may have been concentrated in two

ways. First, the reported intensities are averages over the

entire event so that much higher intensities would prevail

for short time periods. Second, vegetation causes a spatial

redistribution of incoming rain such that small areas of

much higher intensity than the average are expected. On the

other hand, it should be considered that Ksat is the minimum

infiltrability and that lateral movement away from local

high intensity spots is likely.

Rapid changes in water table depth and/or stream

discharge have been cited as evidence for bypassing (Thomas

and Phillips, 1979; Beven and Germann, 1982). Thomas and

Phillips (1979) noted that water in macropores can flow into

or below the rooting zone in a matter of minutes and de-

scribed flow from a spring 30 minutes after cessation of a

large (4.57 cm) rainfall.

From 1 September to 24 July, 1982-1983, the water table

depth was monitored about 5 times each week. At the same

time, a daily record of precipitation was maintained.

Measured water table depths ranged from 80 to 200 cm. The

length of time between rainfall events and rise in water

table depth could be estimated only to within 24 hours

because measurements of both rainfall and water table depth

were made at approximately the same time (7:00 AM). There-

fore, if a rise in water table depth and rainfall were

recorded on the same morning, the water table response could

have occurred between 0 and 24 hours after the rainfall

event. On the other hand, if no change in water table depth

was recorded the following day, either the incoming water

was absorbed in the soil above the water table or the

response took longer than 24 hours.

During the 11 month measurement period there were 27

rainfall events of less than one cm had no measurable effect

on water table depth. Among the rainfall events greater

than 1 cm, there were 27 for which the water table depth was

recorded both the day of the event and the day following and

that were not immediately preceded by large (>1.0 cm)

events. Of these there were 23 for which no response was

recorded the following day. The remaining 4 events occurred

when the soil was relatively moist (there had been between

3.2 and 6.6 cm of rain during the preceding 5 days) and the

water table was between 128 and 143 cm deep. If conditions

of field capacity and moderate (0.5 cm hr-1) rainfall

intensity are assumed, infiltration as calculated by the

Green-Ampt method (see Hillel, 1971) to a depth of 100 cm or

greater within 6 hrs is expected. Given this estimation,

rapid arrival of water at the water table is expected.

Thus, no evidence in support of macropore flow could be


From this information it appears that macropore flow

was not common and that it likely did not occur during the

period of measurement. This is consistent with other find-

ings in the field and laboratory. From this we conclude

that the CD model, or other simplified models based on the

concept that all soil-water participates in convective

transport of solutes, should be appropriate for describing

solute flow in this soil.

It is important to note that, at this time, we do not

have criteria for simple determination of bypassing in the

field. Recent work by Russel and Ewel (1985) performed near

our experimental site reported considerable amounts of flow

through selected channels. Many of the soil characteristics

described above as being related to bypassing on a large

scale were present at that site. These included, steep

slopes, the presence of large coarse fragments, a mixed

canopy, and the presence of native soil fauna. It may also

be noted that their observations of bypassing were restrict

to two large events and relatively small portion of the

total surface area. While the effects may be significant on

a hydrologic scale, they may not be in terms of calculating

nutrient losses by leaching.


Solute breakthrough curves resulting from miscible

displacement of 3H20 in undisturbed soil columns under a

range of soil-water tensions were evaluated in terms of the

mobile-immobile water model and the convective-dispersive

model. Model parameters derived from curve fitting indicat-

ed that the convective-dispersive model accurately described

breakthrough curves performed under tensions greater than

0.1 kPa while the mobile-immobile model better described

breakthrough curves performed under soil-water tensions less

than or equal to 0.1 kPa. Thus, in terms of model selec-

tion, bypassing was significant only when soil-water con-

tents were at or very near saturation.

Dye patterns obtained under saturated conditions showed

that soil-water flow (and thus convective transport) was

confined to small regions within the columns. However, no

easily identifiable, discrete channels were observed in

these regions. It appears that flow was conducted via a

series of relatively large pores, or continuous

pore-sequences. The net effect of the pore sequences on

soil hydraulic conductivity was identical to discrete

channels, but the identification of the channels responsible

is virtually impossible.

Application of tension appears to have disconnected the

most rapidly conducting pore-sequences, thus reducing the

skewedness of breakthrough curves. Under unsaturated

conditions the conducting pore-sequences were slow enough,

and well enough interconnected with the rest of the soil,

that the convective-dispersive model was applicable. Even

so, dye patterns showed that flow was very heterogeneous, as

was also evidenced by high dispersion coefficients.

Comparison of the frequency distributions of

field-measured Ksat values and measured rainfall intensities

indicated that saturated conditions, and hence significant

bypassing, are not expected to occur in this soil.


Observation of water table response to rainfall events

supports this. Based on these experiments, we conclude that

field-scale models based on the convective-dispersive model

for solute movement should be applicable this soil.



There is a large and increasing demand for increased

food production from humid tropical regions. The climate of

these regions is characterized by very high annual rainfall

that is considerably in excess of potential

evapotranspiration. Under these conditions there is a large

potential for loss of nutrients by leaching. This has been

cited as a major limitation to the transfer of fertilizer

technologies developed in temperate zones to humid tropical

regions (Engelstad and Russel, 1975). It is especially

important, therefore, that loss of nutrients by leaching be

considered in the development of cropping systems in the

region. One approach that has been proposed is to manipu-

late cropping systems so that land is continuously planted

to a variety of crops (Cox and Atkins, 1979). Reduction of

leaching loss is among the many benefits ascribed to such


Measurements of nutrient loss from agricultural fields

in the tropics are rare (Greenland, 1977; Keeney, 1982;

Omoti et al., 1983). Those results that have been published

indicate a very wide range in leaching losses (e.g.,

Sanchez, 1976). The lack of leaching studies is due, at

least in part, to the numerous difficulties associated with

such measurements. The fundamental difficulty is that there

is no known means of directly measuring solute leaching

losses. Indirect methods must be employed.

One method that is commonly employed is to use a mass

balance based on the equation:

AF = F u Fl, (3-1)

where F (kg ha-1 ) refers to the change in storage of any

given nutrient and the subscripts "a", "u", and "1", refer

to the amounts of nutrient applied, taken up by the crop,

and lost by leaching, respectively. From an agronomic point

of view, this is reasonable since the amount taken up by the

crop, which is really the object of nutrient application, is

evaluated. The chief limitation to this approach is that

Eq. 3-1 is incomplete. In addition to plant uptake and

leaching, nutrient dynamics are governed by other processes,

e.g., immobilization and denitrification of N and fixation

of K. Another problem is that the amount of nutrient stored

in the root mass is usually poorly known.

Another method used to calculate nutrient leaching is

to calculate the convective solute flux from the the rooting

zone with the equation:

J = qC (3-2)

where J is the solute flux (mg cm-2 day-1), q is the

soil-water flux (cm day-1), and C is the soil solution

concentration (mg cm3 ). This provides a direct estimation

of leaching loss that is independent of other processes that

effect the fate of an applied nutrient. The method requires

knowledge of q and C. The value of q can be calculated

using Darcy's law, but the spatial variability of the soil

hydraulic conductivity is so large that this approach is

generally impractical for field-level research (Warrick and

Nielsen, 1980). An alternative is to calculate q from water

mass balance. This method provides an areal estimate of the

net water movement. Estimation of C can be based on mea-

surements of the soil solution at specific locations in the

field, which can be made with quantifiable precision.

This approach can be extended to trace the movement of

a solute through the crop rooting zone. This provides

information concerning the length of time a solute (e.g., an

applied nutrient) can be expected to reside in the crop

rooting zone under known soil, plant and weather conditions.

Water Balance and Nutrient Leaching

If the system of interest is considered to be delineat-

ed by a unit surface area of soil extending to the depth of

the crop rooting zone, then, assuming that the mass of water

is conserved,AW = W. W where W is the mass of water
in out
and the subscripts "in" and "out" refer to water entering

and leaving the system. This expression can be decomposed

into a number of components such that:

AS + ARW = (P+I+UP+RON) (E+T+APA+ROF+DP), (3-3)

where S = soil profile water content

RW = water held in the roots,

P = precipitation,

I = irrigation,

RON = runon of water from adjacent soil,

UP = upward percolation,

E = evaporation,

T = transpiration,

PA = water stored in the above-ground biomass,

ROF = runoff of water to adjacent soil, and

DP = deep percolation,

with all components expressed in units of cm (the density of
water is assumed to be 1 g cm ). This general equation can

usually be simplified considerably when adapted for a

specific location and set of objectives. Based on observa-

tions in the field, I, RON, and ROF equal 0 in this study.

In addition, the assumption was made that, relative to the

magnitudes of the other components, ARW can be ignored and

APA can be lumped with transpiration. For practical rea-

sons, E and T were lumped into the single term,

evapotranspiration (ET). Thus,

AS = (P+UP) (ET+DP). (3-4)

If the net deep percolation (NDP) is defined as DP-UP and

substituted into Eq. 3-4, the result is as follows:

NDP = P ET AS. (3-5)

In general, S and P can be measured directly, ET can be

estimated from measured parameters, and NDP can then be

calculated as the difference. It is apparent that the

accuracy with which NDP can be determined in this way is

dependent on the accuracy with which the other three parame-

ters are determined.

This fairly straightforward approach is complicated by

the fact that ET is a function of S (or soil-water poten-

tial) when conditions of soil-water stress are encountered,

which requires either that S be determined frequently, or

that S be estimated between measurements. In addition, it

provides no information about water movement between S

measurements or within the system.

Models based on the use of Eqs. 3-2 and 3-5 that

incorporate soil and crop parameters have been developed to

overcome these limitations (Rao et al., 1981; Rose et al.,

1982a). They apply this approach incrementally to layers

within the system and can thereby provide more information

about water and solute movement over depth within the system

and out of the system over time.

Model Development-Concepts

The model used in this study fits the general classifi-

cation of capacity-parameter based models (Addiscott and

Wagenet, 1985). Water movement is calculated from differ-

ences in amounts of water stored. Specific processes

governing the rates of movement are explicitly ignored.

The advantages of this approach over other, rate models, are

that the required parameters are generally easier to measure

and exhibit considerably less spatial variability (Warrick

and Nielsen, 1980). The disadvantage of capacity models is

that, because they are empirical representations of process-

es, they lack the flexibility of more process oriented

models. For this reason it is important to carefully

consider the relationship between the processes of interest

and the empirical representation of them.

There are two basic outputs desired from the model:

(1) NDP, which is to be used to calculate movement of

nutrients beyond the rooting zone using Eq. 3-2 with q =

NDP, and (2) the time of solute residence in the rooting


The following two assumptions are considered fundamen-

tal to the approach: (1) the position of a solute front of

a nonadsorbed, conserved solute is dependent only on convec-

tive transport (i.e., dispersion and diffusion are not

considered), and (2) drainage of soil water proceeds to a

consistent soil-water content called the field capacity

Examination of solutions to the convective-dispersive

equation indicate that, except where v L/D (discussed in the

previous chapter) is very low (i.e., when dispersion is very

high), dispersion has relatively little effect on the

position of a solute front. That is to say that solutes

move in roughly symmetric pulses through the soil. Even

under field conditions, where D is expected to be quite

high, the effect of dispersion on the position of a solute

front is small compared with other errors and the accuracy

that is usually considered acceptable in field studies (Rose

et al., 1982b).

From the first assumption, the traditional

convective-dispersion equation is simplified to

98c/8t = -qac/az, (3-6)

where e is the soil-water content (cm cm ). If one is

interested only in tracking the depth of the solute front

(DSF), the equation becomes;

[3DSF/Dt] = (q/8) = vo, (3-7)

where DSF is expressed in cm, and v0 is the average pore

water velocity (cm day-1). The importance of v0 as dis-

cussed in the previous chapter is apparent. Studies using

packed soil columns have shown that the above principles are

applicable to a wide range of textures and initial condi-

tions (e.g., Dahiya et al., 1984). In general, it is

expected that the approach will be valid for most

field-scale applications if v0 is an appropriate descriptor

of soil water movement. It is clearly inappropriate when

water flows preferentially through large macropores or

around large soil peds or aggregates.

Several studies have demonstrated the validity of the

approach under a wide range of field conditions (Rao et al.,

1976; Cameron and Wild, 1982; Barry et al., 1985). Such

success, however, is by no means universal. Smith et al.

(1984), working with a variety of soils over a three year

period, found that agreement between measured solute front

depths and those calculated with Eq. 3-7 ranged from excel-

lent to poor, varying both with soil and year. Still other

studies have shown the approach to be completely inappropri-

ate (Addiscott et al., 1978; Bouma and Dekker, 1978; Thomas

and Phillips, 1979; White, 1985). In these studies

preferential flow through macropores causing bypassing was

cited as the reason for failure of agreement.

This work demonstrates that the applicability of the

model to the specific conditions of interest should be

tested. In the previous chapter it was demonstrated that,

at the scale of 12.5 cm diameter undisturbed soil columns,

such bypassing did not occur under flow regimes encountered

in the field. Further evidence was cited that such behavior

did not occur at the field level.

The second assumption (concerning fc) has seen a long

history of controversy and it is clear that its application

should depend on the objectives of the study and specifi-

cally measured values. It is based on the general observa-

tion that, upon cessation of infiltration, soils drain with

decreasing rapidity over time. For uniform, well drained

soils, drainage is commonly described with the following

expression: S = at-b, where a and b are arbitrary constants

(Hillel, 1980). Thus, drainage never actually ceases, but

becomes increasingly slow in finite time. The choice of 8fc

and the length of time required to reach that value are

arbitrary and depend on the soil in question and the objec-

tives of the study. In addition, it is important to recog-

nize that the relationship between soil water tension and

soil-water content at field capacity depends on the soil

(Ratliff et al., 1983).

When the second assumption is made, the eventual depth

of leaching from a given infiltration event (I), if ET

during redistribution is negligable, is

DSF = I/fc, (3-8)

where I is the amount of water entering the soil during an

infiltration event (Rao et al., 1976). The significance of

the assumption of a field capacity is that it allows the use

of the amount of the event, which is relatively easy to

measure and use in computations, as opposed to using rates

and time intervals. Equation 3-8 has been modified ro

account for the effects of evapotranspiration and movement

of the solute front within the soil profile (Rao, et al.,

1976; Davidson et al., 1978; Rose et al, 1982b).

Water Table Effects

Application of the modeling approach described above

has been restricted to conditions in which a water table was

not present. The presence of a fluctuating water table in

close proximity to the rooting zone renders the second

assumption invalid because the soil-water does not drain to

field capacity as it is generally defined but to some

greater value that depends on the location of the water

table. In the extreme, at the water table interface, that

value is the saturated water content of the soil (8 ). At

increasing distances above the water table that value

decreases, until at some point, it equals the "normal" field


The soil profile can therefore be divided into two

regions; that above the influence of the water table, and

that within the influence of the water table. The propor-

tion of the soil profile in either region is variable,

depending on the depth of the water table and the height of

water table influence (HWTI). In the region above the HWTI,

drainage is assumed to proceed as if no water table were

present. Drainage in the region below the HWTI is assumed

to proceed to a different, greater value. This value is

considered to be a "temporary" field capacity, 9fc *, that is

dependent on the position of the water table.

The relationship between field capacity, the HWTI, and

the water table depth are illustrated in Fig. 3-1. Note

that the soil-water content at the HWTI is equal to 8fc and

that the water table represents effective saturation. Under

the quasi-equilibrium conditions to which the soil is

assumed to drain, the soil-water tension between the HWTI

and the water table (measured in cm H20) is numerically

equal to the height above the water table (expressed in cm).

This is illustrated as a linear function in Fig. 3-1, but

may be quite different, depending on the soil moisture

characteristic of the soil. This relationship should be

independent of the direction of water table movement if

hysteresis effects are ignored.

6 (cm3cm3)




water table depth

Figure 3-1. Schematic representation of the relationship
between water table depth, the height of water table
influence (HWTI), 6fc*I and 9fc*

Materials and Methods

Cropping Systems and Soils

The two cropping systems chosen for comparison were:

(1) a mixed cropping system composed of laurel (Cordia

alliodora), cacao (Theobroma cacao) and platano (Musa

paradisiaca), and (2) a monocropping system composed of

maize (Zea mays L.). These two cropping systems will be

referred to as the MP, for mixed perennial, and the MA, for

monocropped annual, cropping systems.

The sites of the two cropping systems were separated by

approximately 100 m and were located in the "la Montana"

section of the experimental plots at CATIE (Centro Agricola

de Investigaciones y Ensenanzas) near Turrialba, Costa Rica.

Both cropping systems were planted on Instituto clay loam,

which is classified as a Typic Dystropept, fine, mixed,

isohyperthermic (Aguirre, 1971).

The management level of both systems was designed to

promote high yields and profitability. The MA plot is

located in a larger study area used by Dr. Carlos Burgos to

study the effects of tillage and residue management on maize

yields. It was initiated in November of 1976. A 120-day

variety was planted at a density of 40,000 plants ha-1. Two

crops were planted each year, one in late May, and the other

in early November. Approximately 240, 55, and 40 kg ha1 of

N, P, and K, respectively, were applied each year. Applica-

tions were made on the planting date and approximately 30

days after planting, resulting in 4 applications each year.

The plot dimensions were 32 by 20 m.

The MP plot was part of a larger study of intercropping

with perennial crops directed by Dr. Gustavo Enriquez. It

was initiated in 1977. The planting densities were 1,111;

432; and 123 plants ha-1 for cacao, laurel, and platano,

respectively. The annual fertilizer regime was 140 kg ha-1

N, 30 kg ha-1 P, and 20 kg ha-1 K in four applications. The

dimensions of the MP plot were 18 by 18 m.

Model Developement-Construction

The model used in this study was based on the model,

NITROSIM, developed by Rao et al. (1981). Water movement

was divided into three phases; infiltration, redistribution,

and static. Calculation of soil-water content (9 ) and flux

(q ) for depth increments (i=1,2,...n) of thickness of Az cm

within the soil profile was an iterative process carried out

at discrete time increments (At).

Infiltration was assumed to proceed as a "square pulse"

(i.e. Green-Ampt infiltration) within one time increment at

an infiltration soil-water content ( inf) where 9fc <9 inf
Assumed ET during infiltration was 0. The depth of the

wetting front (dwf) resulting from a rainfall event of I cm

was calculated such that

I = E (einf 9 ) dwf, i=l,2,...,n. (3-9)

The change in DSF (ADSF) resulting an infiltration event, I,

was calculated as

ADSF = 0 I < AW (3-10)

ADSF = (I AW)/8 I > AW (3-11)


AW = E(8inf 8 )Az, i = 1,2,...,m (3-12)

and m is the depth increment in which DSF resided prior to

the event.

Redistribution calculations were based on the following


8 = 9fc + (inf fc) exp(ct), (3-13)

where TRD is the length of time required for the soil to

drain to field capacity. In this way 9 decreases "exponen-

tially" to a value within 1% of 8f after TRD days of


In the algorythm used to update 8, q, and DSF described
i' i
below, 8 is used to denote 8 at t = t + At. Calculations

proceeded from the upper depth increment (i=l) down as


q1 = 0 (3-14)

8 = 8fc + (81 f) exp(cAt) (3-15)

q2= (81' 81) Az/At (3-16)

and, for i > 1 and 9 > 8fc

S = 9 9- 9fc ) exp(cAt), (3-17)
i+l *' i
q = (89 8) Az/At + q, (3-18)

or, for < 9f ,

8 = q Az/At + 9 9 < 98 (3-19)

For depth increments within the crop rooting zone ET is

included so that

1 = 81 U1 At (3-20)

where U1 is the rate of ET from the i-th depth increment.

Changes in DSF (ADSF) were calculated as

ADSF = DSF + q1 (At/Aze ) (3-21)

when DSF was in the i + 1 depth increment at time t.

Some modifications were required to incorporate water

table effects. The water table depth was an input to the

model taken from measured and interpolated values. The HWTI

and 9fc (Fig. 1) were calculated on a daily basis from the

water table depth. All water entering depth increments

within the HWTI in excess of that required to account for

measured changes in water table depth was considered to be

NDP. Redistribution within the HWTI proceeded as described

above with 8fc replaced by 8fc. Upward movement of water

was not explicitly considered.

The static phase commenced after field capacity was

achieved. During that phase changes in soil-water content

occurred only as the result of ET. Note that redistribution

could result from either infiltration events of a lowering

of the water table.

The soil profile was divided into 22, 5 cm increments

for numerical computations in the simulation model. Thus

the total depth of interest of 110 cm. The time step used

for calculation of water and solute movement was 0.1 day.

Days were considered to begin at 7:00 AM, when measurements

of water table depth, rainfall, and soil-water content were

made. Rainfall events were assumed to take place at 4:00

PM, which is roughly the most common time for such events.

All programming was done on a VAX-11 minicomputer in


Determination of Inputs

Use of the model requires as input the rainfall,

potential evapotranspiration, water table depth, initial

soil water content, and the soil and plant parameters.

Rainfall was measured daily with gauges located adjacent to

each plot. The water table depth to 220 cm was measured

about 5 times each week in perforated PVC tubes located in

the center of each plot that served as piezometers. Model

outputs were generated for two different time periods, the

simulation period and the calibration/validation period,

which are described subsequently. The initial soil-water

content for the simulation periods was determined

gravimetrically and for the calibration/validation period

with a neutron probe.

Soil samples for soil-water content determination were

taken at 10 cm intervals from 6 auger holes in the MP plot

and 9 holes in the MA plot. Soil-water content was deter-

mined by weight loss after drying at 1100 C to a constant

weight. This was converted to a volumetric basis by multi-

plication with the soil bulk density.

Two methods were used to calibrate the neutron probe.

First, counts were made simultaneously with gravimetric

sampling approximately 1 meter from the access tube.

Second, a large tub was filled with soil taken from an

adjacent plot. The soil was then mixed and adjusted to

three different soil-water contents, packed into the tub,

and measured. Counts measured in the tub were adjusted for

bulk density as suggested by Greacen et al. (1981). The

resultant calibration curve, including field and

tub-measured values are shown in Fig. 3-2. There were five

access holes in each plot and counts were made at 15-cm

intervals from depths of 20 to 110 cm.

Evaporation from a class A US Weather Bureau pan

located approximately one km from the site was the basic

input used for calculating the potential evapotranspiration

(PET). The measured pan evaporation (E pan) was converted to

PET using the equation PET = k Epan where k is the pan con-

stant, taken from Doorenbos and Pruitt (1974) as 0.8.

Soil Parameter Estimation

Soil profiles were assumed to be uniform over the

depths of interest with respect to the required soil hydrau-

lic parameters because only small trends with depth were

observed in the plots. There were some differences in soil

properties between the two plots used in the study even

though they were separated by only 100 m and mapped as the

same soil series. The main difference was that soils in the

MP plot had about 5% more clay, on the average, than the MA

plot, which corresponded with noticeably higher soil-water

contents. For this reason slightly different parameter

values were used for the two plots.


00 0











o Tub calibration

* Field calibration

1.3 Z




0.30 0.35 0.40 0.45 0.50 0.55 0.60


Figure 3-2. The relationship between neutron probe count
ratio and soil-water content used to calibrate the probe.
The regression equation was; count ratio = 2.69 9 + 0.65,
with r = 0.993.

.0 *

The soil parameters required to calculate infiltration,

drainage, and movement of the solute front in depths above

the HWTI are 8 inf 8fc, and TRD. Sixteen tensiometers,

placed at 15 cm intervals from 10 to 115 cm deep, 2

tensiometers at each depth, were installed in a 3x3 m

subplot immediately adjacent to the MA plot. Water was

ponded until tensiometers indicated that approximate steady-

state flow conditions had been achieved. At that point the

subplot was covered with straw mulch and plastic to elimi-

nate evaporation, and changes in soil-water content were

monitored with a neutron probe. The value of TRD was

estimated as the length of time required for the rate of

change of soil-water content to become negligible.

A similar test was attempted in the MP plot but had to

be discontinued after large rainfall events caused a rise in

water table. The value of 8fc was estimated as the soil

water content at the soil-water tension (cm) equal to the

HWTI (also in cm). The soil water characteristic curve

under desaturating conditions in two 5.4 cm diameter cores

of undisturbed soil from the MP plot was measured for this

purpose. The value of TRD was assumed to be the same for

both plots. Estimation of 8 was based on observation of

the soil water content over time in both plots.

Two additional soil parameters are required for calcu-

lations within the zone of influence of the water table; the

height of water table influence (HWTI), and 9 fc Average

measured soil-water tensions at two different depths, 75 and

105 cm, were compared with the water table depth. The

soil-water tensions were measured with tensiometers at the

two plots. In the MP plot there were 10 tensiometers at 75

cm and 5 at 105 cm. There were 9 tensiometers at both

depths in the MA plot. The HWTI was taken as the maximum

elevation difference for which approximate equilibrium

between soil-water tension and the water table could be
expected. Values of 9fc were assumed to vary linearly with

the height above the water table as shown in Fig. 1.

Plant Parameter Estimation

Measurements of root weights in the MP system shown in

Table 3-1 (Alpizar, personal communication) were used to fit

the following exponential relationship of root concentration

with depth:

Root weight = 25.179 exp (-0.067z), (3-22)

where z is depth (cm). The maximum rooting depth was taken

as 70 cm. This rooting depth was assumed to be constant

over time.
Table 3-1. Root distribution in the MP plot .

Depth Mean Root-Weight Standard Deviation
cm g cm
0-15 14.0 8.7
15-30 6.4 3.3
30-45 0.4 0.4
*Alpizar (Personal Communication)

No measurements of root concentrations were made in the

maize plot. Root distribution with depth over time was

calculated from relationships observed by NaNagara et al.

(1976) as applied by Davidson et al. (1978). The rooting

depth during fallow periods was assumed to be 35 cm.

Estimated PET values were converted to estimates of

actual evapotranspiration (AET) using a cropping factor (CF)

in the following equation: AET = (CF)(PET). A constant CF

value of 1 was used for the MP plot. This consistent with

estimates for cacao grown alone (Doorenbos and Pruitt, 1974)

and the fact that canopy coverage on the plot was complete.

The CF value for the MA plot varied over time as the

maize crop was planted and matured. Estimates of CF during

the time that the maize crop was present were taken from

Doorenbos and Pruitt (1974), who divided the cropping season

into 4 stages with respect to water use. Estimation of CF

in the MA plot during fallow periods, which had considerable

weed growth, was made by calibration as will be described in

the next section.

The amount of water extracted from each layer was

calculated using the approach of Molz and Remson (1970).

Estimation of soil-water stress was based on the following

relationship which is similar to that used by Davidson et

al., (1978):

AET = AET (8-8 )/(8 -e ), 9 < 89stress' (3-23)
cal pwp stress pwp stress (
where AETcal is the transpiration calculated in the absence

of water stress and stress and 8 are the soil water

contents of stress initiation and permanent wilting, respec-



Soil-water measurements were made periodically with the

neutron probe during the calibration/validation period.

There were 5 neutron probe access tubes in each plot. The

model output used for these purposes was the profile water

content, S. This is defined as

S = fedz (3-24)

where z is depth and the lower boundary is at 110 cm. This

was calculated from measured 8 values as

S = E9.Az, n=1,7 (3-25)

where i represents different depth increments, and Az the

increment thickness. As a check on the averaging procedure

used, S values determined by gravimetric sampling and

neutron probe measurements 18 July were compared.

For logistic reasons, the calibration/validation period

was initiated at different dates at the two plots. It

started on 14 April in the MA plot and 18 April in the MP

plot. These data were used for two purposes, calibration

and validation.

The first 18 days of the calibration/validation period

were used for model calibration. There was no rain during

that time, which was preceded by 3 weeks in which only 4.6

cm of rain were recorded. Thus, DP and P during those 18

days were 0.

Soil-water tensions of greater than 60 kPa by 23 April

and changes in S that were less than that expected indicated

that the crops on the MP plot probably experienced water

stress during the period. The parameters stress and 8 p

were adjusted to account for this, assuming that changes in

S were due to AET. Soil-water tensions on the MA plot were

much lower at that time (about 22 kPa) due to the reduced

transpiration of the fallow vegetation. Measured changes in

S during this time were used to determine the cropping

factor for the fallow vegetation, again assuming that all

changes in S were the result of AET. The water stress

parameters for the MA plot were taken from characterization

data of Aguirre (1971) but were not critical because dry

periods did not occur while the maize crop was actively


After 2 May there was consistent rain and the condi-

tions of P and DP equal to 0 did not hold. From this point

on the calibration ended and measured S values were used to

evaluate the accuracy of model calculations.


Simulation of water and solute front movement began on

4 Nov. 1982 for the MA plot and 24 Nov., 1982 for the MP

plot and ended for both plots on 18 July 1983. Initial

soil-water contents were determined gravimetrically. The

starting date for the MA plot corresponded to the maize

planting date. A second planting of maize in the MA plot

took place on 30 May, 1983. All parameter values used in

the simulation were determined as described above.

Net Leaching Loss

Net leaching losses over the simulation period were

calculated using Eq. 3-2 with the NDP used in place of q and

C representing measured soil solution concentrations. Soil

solution concentrations of NO3 NH 4 Ca2+, Mg2+, and K

were measured in extractions from 7.6 cm diameter ceramic

porous cup samplers located at a depth of 90 cm. The cups

were pretreated in 0.1N HCL prior to installation. There

were 9 samplers in the MA plot and 8 in the MP plot.

Collections were made at approximately biweekly intervals

starting 30 Nov. in the MA plot and 11 Dec. in the MP plot.

Samples were collected by applying 40 kPa tension overnight.

No collection was made when the soil-water tension was

greater than 40 kPa. Concentrations of NO3 and NH4 were

measured using a steam distillation technique (Bremner,

1965) and the cations were measured by atomic absorption.

Depth of Solute Front

The depth of the solute front movement resulting from

application of a nonadsorbed nutrient (e.g., NO3 ) was

simulated for 3 application dates in the MA field plot and 2

in the MP plot. The 3 application dates chosen for the MA

plot, 2 Nov., 30 Nov., and 30 May, were actual fertilization

dates on that plot. Calculations for the MP plot at the two

latter dates were included for purposes of comparison.


Model Inputs

The depth to the water table and measured rainfall

amounts over the simulation period are shown in Fig. 3-3.

The water table was consistently deeper under the MA plot

than under the MP plot. A slight elevation difference

between the two plots (the observed water table differences

E 8.0

* 4.0


OSV u 'MPplot

a 14O



230 1 I .
30 60 90 120 150 180 210 240
Time (day)

Figure 3-3. Precipitation and water table depth at both
plots over the simulation time period (2 Nov., 1982 to 18
July, 1983).

were about 20 cm) and differences in the drainage system at

the two plots, explain this. Open ditch drains of 100 to

150 cm depth drained both plots, but the drains around the

MA plot had a greater gradient and removed water more

quickly to a greater depth. After that depth was achieved,

both plots behaved similarly.

The rainfall distribution was fairly uniform during the

early portion of the simulation period and was then charac-

terized by two relatively severe droughts. Shortly after

day 180 rainfall was consistently high.

Measured pan evaporation varied slightly during the

simulation period. It ranged from 0.26 cm day-1 in Dec. to

0.48 cm day-1 in April. The average monthly pan evaporation

over the simulation period was 0.35 cm day-1.

Parameter Estimation

A summary of parameter values obtained is shown in

Table 3-2. The change in soil-water tension with depth over

time during the determination of 98f and TRD is shown in

Fig. 3-4. It indicates that the profile is roughly uniform

with depth with respect to soil hydraulic properties. The

corresponding changes in S are shown in Fig. 3-5. The

overall course of drainage is similar to that reported

elsewhere in the literature (Hillel, 1980) in that In S is

proportional to In t. This relationship was approximated

using the following parameter values: 9inf = 0.49, c =

0.46, and TRD = 5.0 days. The curve calculated using the

above parameters is also shown in Fig. 3-5. The calculated


z -2.9


S4.0-0 ------

20 40 60 80 100 120

DEPTH (cm)

Figure 3-4. Soil-water tension as a function of depth and
time after ponding. Time after ponding, in days, is
indicated below points representing the upper depth.


54.0 fitted (smooth curve)

53.0 *



50.0 I
I 2 3 4 5 6 7 8

TIME (days)

Figure 3-5. Profile soil-water content after cessation of
ponding. Fitted values were calculated with model
parameters in Table 3-2.

values are close to those measured for times up to 7 days.

This is considered to be sufficient because drainage rates

for longer time periods are expected to be considerably

different due to changes in gradient that result from plant

uptake of water and evaporation.

Table 3-2. Summary of model parameter values.
Cropping System
Parameter MA MP

S0.46 (cm tcm-3) 0.52 (cm cm-3)
9-i 0.49 0.55 "
9 s 0.55 0.57 "
9sat 0.29 0.38 "
8pwp 0.42 0.48
THress 5 (day) 5 (day)
Max. Root 100 (cm) 70 (cm)
CF variable 1.0

The measured soil water characteristic curves upon

which the estimation of 9fc in the MP plot was based are

shown in Fig. 3-6. The value corresponding to a soil-water

tension of 5.5 kPa (or 55 cm of H20), 0.52, was chosen. The

tension of 55 cm H20 is numerically equal to the value of

HWTI chosen (see discussion below). The estimated values of

8 and inf are shown in Table 3-2. The TRD was assumed

to be equal for both plots.

The average measured soil-water tension and distance

from the tensiometers to the water table in both plots are

plotted in Figs. 3-7 to 3-10. These plots were used to

estimate the HWTI. Each plot is somewhat distinctive, which

is due, in part, to the differing conditions encountered at

the different depths in the two plots. Recalling Fig. 3-1,

the HWTI should be the maximum distance between the water



0 a-)

-1 0

00 1
r 0

4 -r-1
O <
0 /"1 .0
-- t-4

M E-

| /" (O : :
)) 0O

0 0 a


/ >H -
o 0

D +4J U)



IC) I- L- <
(0 -1 >

O Ou I) In

( w W1O) 6D 43)


* S

**0 0
* *
* *
* *5

\ A.
\ *



. Ir
0 x



0 >



(Pd A) NOISN31 39VJ3AV



H0-I 0 ()
0) 0

0 *

-Hr r
- 0Z
4-J Z g
o a

0 CO (a
4-)-' +
4 ;aOn C O
a0 -i-

I-p4Q) C
0) (1) 0)
4+) 4J r4-
0 aU

+ 4J 0)(
0H ) 0
a a)

0O 0P -4
S Q)co
*H 0 U

0) Cz
4 4.1-1 01

*H *04 -1
0) 0 iH *c

JQ 0)4

r Irm

* *
0 0

* *0

\ '

\ ** *





O _


o 5








aa )

02 4) (4-
r M a)
4-) -H
0 0z

(U 4Q
m a-H
4- 04

01 4J

(u c

ri a) -ra)
-ri C a a)
0 -
w 0 r

-P U C5

(0 Clu 4-
a) aW --
O U -H
0 412

a CDQ 4-
a*0 aiQ)
0 r-C CT

-4 (a)

C4 0 -
Q) 4J 0

4 0Q) H r-I

OH 4
4- U W0 4

c) 41 4 (u

)C *4

U rIh *
(l O 4
n S, Pl

* *

(vd >1) NOISN31i 39V3AV



0 0e
-,-D 4

4 rO


_* 0o (1) i
%* 0 *a O r 4
*o W2 0 *
Z 0* 4J .U.

0 .0 O .* o- u-I C0
O co W
0g0 LLJ *3 I iW )
CC > c 0

0 r. 0

.a)C 0U
\ ro w (1

) >\ O ZC

P- 4 rd
qq- / 0 0-H

U) aO 4 U
S- -- u
*C L4 0
Yra 41 a) -r
4- rd4 rd
ooo oEn -P 4-)

\( d >1 S sV A
a ) r;) 04

\O | LO I0 O
N C-j c 4r 43

\ (0 4 r }-4W
(Bd A) NOISN31 39)Va3AV !| H

4 Q) U)

a) F.
r- 3 -P i

E 0 r.

0 a)
\* O C

a) 0 ) 0

M ** 0 w
\ E* C) t

** 0 > lj- 4J4J
\ ) 4-J 4j 4

\7** z 0 0 *
'- C (uT C


4 s-!J 3
\ a in a) s
S- C)o
I- C) ""I

0 0 0
*:\ r: 4 4-

** :4-) C)- E

a) Ha)
S\ 0 ()J U -r

\;I L-N r.

6 0\0 0 o o<
S4 .Q 0a)

table and tensiometers at which equilibrium can be assumed.

When that distance is greater than the HWTI, soil-water

tension could either be greater than equilibrium due to ET,

or it could be less, due to the attainment of field capacity

or infiltration. Of course, equilibrium values could also

be observed.

Considering the errors that could be introduced by

varying elevations within the plots and the accuracy with

which measurements were made, reasonable estimates of HWTI

range from 45 to 65 cm. A value of 55 cm was chosen for

HWTI for modeling purposes, which should be approximately

correct in all cases. This value agrees reasonably well

with the tension measured at the plots after 5 days drainage

(Fig. 3-4). Discrepancies on the order of 1 to 2 kPa are

not considered important because they have little effect on

estimated 9 values (Fig. 3-5). Although the tension at

field capacity is somewhat lower than commonly assumed, it

is similar to values measured in other studies (Wierenga,



As a check on the overall averaging procedures used to

calculate S in the field from measured Q values, S calculat-

ed from 0 measured gravimetrically was compared to S calcu-

lated from 9 measured with the neutron probe. The results

shown in Table 3-3 indicate good agreement.

Table 3-3. Comparison of profile water contents measured
gravimetrically and with the neutron probe.
MA Plot MP Plot
Method Ave.S S.D. C.V. n Ave.S S.D. C.V. n
cm cm
Grav. 50.99 2.30 0.045 9 55.48 2.59 0.047 5
Probe 50.79 0.73 0.014 4 56.62 1.76 0.031 4

The value of CF on the MA plot during the fallow period

obtained by calibration was 0.5. This value was used in the

remainder of the validation run and in the subsequent

simulation run for fallow periods. The other values ob-

tained by calibration were 8 stress and 8 P for the MP plot

(Table 3-2).

Measured values of S over time in both plots are shown

in Figs. 3-11 and 3-12 along with simulated results. There

was close agreement between simulated and measured values on

both plots. The variability of S estimation as indicated by

the 80 and 90% confidence intervals (vertical bars in Figs.

3-11 and 3-12) in the MA and MP plots, respectively, was

rather high. This was not due to high variability per se

(the coefficient of variation ranged from 1 to 5%), but to

the limited sample numbers. Efforts to improve that condi-

tion were hampered by defective aluminum tubing.

Soil Solution Concentrations

Measured soil solution concentrations over time under

the two cropping systems are shown in Figs. 3-13 and 3-14.

In general, the concentrations under the MA system are at

the lower end of the range reported by Barber (1984) for

concentrations in the surface horizons of Midwestern soils.

By that standard, the concentrations in the MP plot were



QJ r-H

0 oo0

o .4 f
OR 0

o co

a= E 0 .o
-/ 2>cn ^O >
S t() *-

iH -H
,5 J ,, O

o con

e or*

r *,---H <

o ^ o


o 0 rd


-4- t o *

-./ r- ( q

7- U) > 0


o 0 00
-4O 4
c0 o 04

a OD 0

4a- 0 0

0 zJ

S> (0

0 0Q

o c-
o 0E
C o-o

-4 -0 0
N0 >
o H > o

f to O N LOo 0 ? 4 -r

n 0 *H4
S ) S > w *

-4 0E

'I I u




co q** 0 CD 0 d* 0
0 0 0

(NW) N


+ C
o -


- .0
+ 4J


m0 H
41 0


(O *r4

U -H -

*ri -

l *44

4^ a rdO
3 3 <
3tji 1
**- *
fa( Mr 4^

(VYW) )i

r-1 01\



M -H

- 1-I


.,- (0



0 04

O U,


U 0


UI r-o
,( 0
r-l c l


C 0

P4 4-

0 0 0 O 0

0 0 6 O0

(vW) 6qW

(V~w) .80

quite low for all ions except K Concentrations of

NH4 were relatively low under both systems and were not

included because the concentrations measured were generally

at or below the minimum concentration required for accurate

determination for the methods of analysis used.

With the exception of K soil solution concentrations

under the MA plot were clearly greater than those under the

MP plot. Even in the case of K the soil solution concen-

trations in the MP plot were consistently lower than those

in the MA plot.

Examination of the concentrations and the associated

90% confidence intervals, suggests that the concentrations

were fairly constant over time. Analysis of variance

indicated that there was no significant difference

(a=0.9-0.99) between the mean concentrations measured on

each plot for all ions measured except NO3 on the MP plot.

This indicates that the means can be pooled. The overall

means and pooled standard deviations are shown in Table 3-4.

Comparison of the means of the ions in the different plots

were significantly different (a=0.99) for K Ca 2+, and
Mg The concentrations of NO3 in the MP plot were more

than one order of magnitude less than in the MA and fre-

quently below the accurate detection level of 0.011 mM.

Table 3-4. Average soil solution concentrations.

Cropping System
Nutrient MA MP
ion Ave. S.D. Ave. S.D.
NY3 0.62 0.349 *
K 0.11 0.074 0.043 0.028
Ca2+ 0.16 0.104 0.016 0.010
Mg 0.13 0.057 0.024 0.006
*Could not be pooled.


Values of S calculated for the simulation period are

also shown in Figs. 3-11 and 3-12. After approximately 150

days of simulation, calculated values agreed well with those

measured. The relatively low S values calculated for the

simulation run were probably due to one of two causes:

either the AET in the dry period prior to the calibra-

tion/validation period was overestimated due to the fact

that the laurel trees shed their leaves during the dry

season, or there was unaccounted upward movement of water

prior to that time.

The greater AET from the MP plot (3-5) was due to the

fact that transpiration from that plot was continuous. This

resulted in simulated differences in NDP (Table 3-5).

Net Leaching Loss

The net leaching loss (NLL) of the four elements

measured were calculated by multiplying the NDP and the

average soil solution concentration. The results are shown

in Table 3-5. These results are for a time period somewhat

less than one year, which makes direct comparison with most

other measured results difficult. The range of values

reported, however, is so great that these values easily fall

within it (Sanchez, 1976). From an agronomic viewpoint, the

loss of NO3 -N in the MA plot is clearly important.

Table 3-5. NDP and Net Leaching Loss Values.

Value Cropping System
AET 45 (cm) 55 (cm)
TR* 111 (cm) 111 (cm)
NDP 66 (cm) 57 (cm)
NLL N 57 (kg/ha) 1 (kg/ha)
K 2+ 3 (kg/ha) 1 (kg/ha)
Mg2+ 21 (kg/ha) 3 (kg/ha)
Ca 43 (kg/ha) 3 (kg/ha)

*TR is the total rainfall during the time of


Depth of the Solute Front

The simulated movement of a nonadsorbed solute front is

shown in Fig. 3-15. The rate of movement from the third

application is considerably greater than that from the

first. Movement is clearly related to patterns of rainfall.

Note that the DSF from the first two applications tend

to converge with time. Also note that little residual

effect of fertilizer can be expected between the two appli-

cations as the front has moved well beyond 110 cm depth

before the second planting. The difference between the MP

and MA plot is ascribed to the greater AET from the MP plot.


c 40- N
LL ----------

U) MAplot -
SMP plot --- ot
o 80-


0 30 60 90 120 150 180 210 240 270

Time (days)

Figure 3-15. Depth of solute front movement resulting from
application of a nonadsorbed solute on days 1, 30, and 207
in the MA plot and days 30 and 201 in the MP plot.

.1 1 0. Ah .111h