Title Page
 Table of Contents
 List of Tables
 List of Figures
 Review of literature
 Methods and materials
 Results and discussion
 Summary and conclusions
 Future work
 Nitrate and iodide analysis using...
 Evaluation of tension lysimete...
 Coefficients for redistribution...
 Literature Cited
 Biographical sketch

Title: Spatial variability of solute leaching in a sandy soil
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00090191/00001
 Material Information
Title: Spatial variability of solute leaching in a sandy soil
Series Title: Spatial variability of solute leaching in a sandy soil
Physical Description: Book
Creator: Flaig, Eric George,
 Record Information
Bibliographic ID: UF00090191
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001638773
oclc - 24160949

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








I would like to express my appreciation to Drs. Ken Campbell and Del

Bottcher for their continued support through the long duration of this

project. Particular appreciation is extended to Dr. Campbell for

editorial assistance during the final preparation of this document.

Financial support from Dr. Bottcher for this study is particularly

appreciated. The financial support from Dr. Mansell during the first year

of my program is appreciated. I would like to thank Dr. Suresh Rao who

introduced me to the concepts of spatial variability and who provided many

probing questions during the early phase of this project. I would also

like to thank the other members of my supervisory committee, Dr. W. C.

Huber and Dr. J. W. Jones. Appreciation is also extended to Dr. W. Graham

for participation in my defense.

Special recognition is reserved for my wife, Beatrice, without whom

completion of this project would have been impossible. The considerable

support she provided in the field, laboratory, and home were essential to

complete copious analysis required to complete this work.

Finally, support provided by South Florida Water Management District

for final completion of the dissertation is recognized.



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

LIST OF TABLES ...................... ............................... v

LIST OF FIGURES .................................................... viii

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

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

REVIEW OF LITERATURE................................................. 5

Introduction.............................. ....................... 5
Solute Transport.. .............................. ................. 5
Solute Transport Modelling ........................................ 12
Soil Water Flux .............. ........... ......................... 23
Sampling Methods..................................................... 30
Spatial Variability of Soil Properties .......................... 35
Spatial Structure of Soil Properties............................. 37
Spatial Analysis ..................... ... ......................... 41
Uncertainty ............................ ......................... 59
Sampling Stategies......................... ... ..................... 62
Synthesis .............................. ......................... 67

METHODS AND MATERIALS ............................................... 70

Introduction.. ........................... ........................ 70
Laboratory ............................. ......................... 72
Field Experimentation............................................. 73
Solute Leaching Loss............................................ 86
Spatial Analysis ................... ......... ......... ............ 97
Uncertainty ....................................... .............. 102
Sampling Strategies............................................... 103

RESULTS AND DISCUSSION .............................................. 106

Application Uniformity ............................................ 106
Solute Concentration Time Series................................. 108
Solute Pulse Velocity ............................................. 117
Solute Mass Loss ................................................ 124
Solute Mass Recovery ........ ...... ................................ 128

Spatial Structure ............................................... 129
Variance Partitioning ............................................ 149
Uncertainty. ............................ .......................... 151
Areal Estimation ................................................. 155
Sampling Strategies ............................................... 157
Application. ............................ .......................... 168

SUMMARY AND CONCLUSIONS ............................................ 172

Summary .......................................................... 172
Conclusions. ............................ .......................... 173

FUTURE WORK ......................................................... 176

Improved Description of Solute Transport......................... 176
Impact of Cultural Practices on Solute Transport................. 177
Sampling Patterns for Estimating Short-Range and
Long-Range Variability ......................................... 177




LITERATURE CITED .................................................... 205

BIOGRAPHICAL SKETCH ................................................. 225


Table Pag

1 Average value (Ave) and standard error (SE) for
initial water content (ei), and residual water
content (8,), bulk density for field plot..................... 75

2 Soil physical properties of typical Kendrick sand
(Loamy, siliceous, hyperthermic, Arenic, Paleudult)........... 75

3 Irrigation and chemigation parameters for field experiments... 81

4 Sample collection times for each experiment in hours past
cessation of chemigation ...................................... 84

5 Location on experimental site (X), Equivalent depth of
irrigation (Depth), Solute pulse velocity (SPV)............... 95

6 Sampling patterns used to evaluate spatial structure of
simulated two-dimension regionalized variable. ................ 104

7 Summary statistics for application: equivalent depth of
irrigation (Depth), total solute applied (Tots), and
concentration of chemigation solution (Co) for Experiments I,
short transect; II, grid; III, transect....................... 107

8 Statistical properties for estimated solute pulse velocity
(SPV, cm hr-1) for three field Experiments: I, short
transect; II, grid; III, transect............................. 118

9 Statistical properties for estimated solute mass loss (SML,
kg ha-1) for three field Experiments: I, short transect; II,
grid; III, transect .................. .................... ..... 127

10 Summary statistics for fraction of applied solute mass
recovered (FSR) for Experiments I, short transect; II,
grid; III, transect ........................................... 128

11 Coefficients of variation (CV) for application variables:
equivalent depth of irrigation (Depth), Total solute applied
(Tots), concentration of chemigation solution (Co), and
solute mass loss variables: solute pulse velocity (SPV),
solute mass loss (SML), and solute fraction recovered (FSR)
for Experiments I, II, and III................................. 131

12 Regression relationships of application: depth of irrigation
(Depth), total solute applied (Tots), solute concentration in
chemigation (Co), and solute mass movement variables: solute
pulse velocity (SPV), solute mass loss (SML), and solute
fraction recovered (SFR) against spatial coordinates for grid
experiment and transect experiments........................... 133

13 Regression relationships used to account for drift for solute
pulse velocity (SPV) and solute mass loss (SML) against
extrinsic sources, depth of irrigation (Depth), total solute
applied (Tots) for SML regression residuals (SML-res) against
intrinsic source SPV for Experiment II, grid, and Experiment
III, transect................................................. 137

14 Semivariogram models for SPV and SML adjusted for applied
irrigation and solute mass for Experiment II, grid and
Experiment III, transect ..................................... 138

15 Partitioning of variance among regression relationships and
semivariogram models for Experiment II, Grid and Experiment
III, Transect................................................ 150

16 Uncertainty in solute pulse velocity (SPV) and solute mass
loss (SML) due to parameter variability for Experiment II.... 154

17 Areal estimates and standard error for solute pulse velocity,
SPV, and solute mass loss, SML, for Experiment II, grid, and
Experiment III, transect calculated by simple average and
block kriging ................................................ 156

18 Confidence limits associated with identification of the
correct semivariogram model from experimental semivariograms.
Confidence limits are presented for lag 4 results of 50
subsamples for simulated data for stationary covariance,
covariance with a random nugget variance, and covariance
with a nugget variance and trend.............................. 160

19 Recovery of iodide and nitrate from porous ceramic cups
soaked in solution of 5-, 10-, and 20-mg 1-1 KNO3 and Nal
made up in tap water and deionized water...................... 184

20 Iodide concentration in consecutive deionized water washes
of porous ceramic cups following uptake of 20 ml of iodide
solution as affected by aging................................. 186

21 Average percent recovery for iodide and nitrate for initial
sampling of 5-, 10-, and 20-mg/l solutions of KNO3 and Nal
in tap water and first deionized water wash of porous
ceramic cups following uptake of 20 ml of iodide solution
as affected by aging .......................................... 186

22 Coefficient of determination, r2; estimated pore water
velocity, v; from CXFIT, and the coefficents, a and b,
from estimation of soil water flux for Experiment I........... 199

23 Coefficient of determination, r2; estimated pore water
velocity, v; from CXFIT, and the coefficents, a and b,
from estimation of soil water flux for Experiment II........... 199

24 Coefficient of determination, r2; estimated pore water
velocity, v; from CXFIT, and the coefficents, a and b,
from estimation of soil water flux for Experiment III.......... 202



Figure Page

1 Soil water during redistribution for Eustis sand following
ponding: a) soil water flux, b) soil water content............ 29

2 Theoretical semivariogram: a) components, b) semivariogram
models. ................................. ....................... 40

3 Experimental site layout. a) Areal view of lysimeter placement
for all experiments, b) Areal view of lysimeter placement for
Experiment I, c) Areal of lysimeter placement for Experiment
II, d) Areal of lysimeter placement for Experiment III........ 78

4 Sampling manifold and collection tubes for tension lysimeters. 83

5 Soil solution nitrate concentration time series and fitted
model from CXFIT with coefficient of determination (r2) for
model fit for three lysimeters (identified by E-W coordinate
and N-S coordinate) for Experiment I: a) Location 14.25, 2:
r2 = 0.96, b) Location 15.50, 2: r2 0.70, c) Location
13.50, 2: r2 0.35............................................ 89

6 Soil solution iodide concentration time series and fitted
model from CXFIT with coefficient of determination (r2) for
model fit for three lysimeters (identified by E-W coordinate
and N-S coordinate) for Experiment II: a) Location 3,5:
r2 = 0.99, b) Location 2,7: r2 0.72........................ 90

7 Initial and final estimated soil water flux relationships
during redistribution for location 18.00, 2 for Experiment I.. 96

8 Spatial distribution of applied nitrate and irrigation during
Experiment I: a) areal application rate of nitrate, b) total
equivalent depth of irrigation............................... 109

9 Spatial distribution of applied iodide and irrigation during
Experiment II: a) areal application rate of iodide, b) total
equivalent depth of irrigation............................... 110


10 Spatial distribution of applied nitrate and irrigation during
Experiment III: a) areal application rate of nitrate, b) total
equivalent depth of irrigation............................... 111

11 Nitrate and iodide concentrations during infiltration and
redistribution for Experiment B8: a) Site B8-5, b) site B8-2.. 113

12 Distribution of solute concentrations in soil solution
collected at z 50 cm depth with time during each experiment;
median, upper and lower quartiles.

sites samples
a) Experiment I 20 17
b) Experiment II 98 12
c) Experiment III 96 20......... 115

13 Spatial distribution of estimated nitrate leaching for
Experiment I: a) nitrate pulse velocity (SPV), b) nitrate
mass loss (SML) ............................................... 120

14 Spatial distribution of estimated iodide leaching for
Experiment II: a) iodide pulse velocity (SPV), b) iodide
mass loss (SML)............................................... 121

15 Spatial distribution of estimated nitrate leaching for
Experiment III: a) nitrate pulse velocity (SPV), b) nitrate
mass loss (SML)................................................ 122

16 Semivariogram for solute pulse velocity for Experiment II:
a) experimental semivariogram, b) experimental semivariogram
and fitted model .............................................. 135

17 Semivariogram for solute pulse velocity for Experiment III:
a) experimental semivariogram, b) experimental semivariogram
and fitted model............................................... 136

18 Experimental semivariogram for solute mass loss (SML) for
iodide in Experiment II ....................................... 140

19 Experimental semivariogram for solute mass loss (SML) for
nitrate in Experiment III..................................... 141

20 Experimental semivariogram for solute mass loss residuals
(SML-resid) adjusted variability in irrigation for iodide
in Experiment II: a) Experimental semivariogram,
b) Experimental semivariogram and fitted model ............... 142

21 Experimental semivariogram for solute mass loss residuals
adjusted for variability in irrigation depth for nitrate
in Experiment III ............................................. 144

22 Experimental semivariogram for solute mass loss residuals
(SML-resid) adjusted variability in solute pulse velocity
for nitrate in Experiment III: a) Experimental semivariogram,
b) Experimental semivariogram and fitted model ............... 145

23 Semivariograms for solute leaching considering different lag
distances: a) solute mass loss (SML), b) solute pulse velocity
(SPV ).......................................................... 147

24 Semivariograms for mean and 95% confidence limits for 50
sample sets of increasing size taken from a simulated
2-dimensional random field. The random field has a known
Bessel covariance with range 8, sill .24, nugget 0.
The simulated random field is 100 x 100; 10,000 points.
Experiment semivariograms are plotted with analytical model:
a) n 5 x 5, b) n 10 x 10, c) n = 20 x 20.................. 162

25 Average error for 50 sample sets of increasing sample size for
different patterns taken from a 2-dimensional random field:
a) mean average-error, b) variance of average-error........... 164

26 Mean squared-error for 50 sample sets of increasing sample
size for different patterns taken from a 2-dimensional random
field: a) mean mean squared-error, b) variance of mean
squared-error................................................. 165

27 Volume of soil water extracted for varying soil water content
for Eustis sand hand packed in a bucket for bulk density =
1.51 g/cm3 .................................................... 191

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



Eric George Flaig

December, 1990

Chairman: Kenneth L. Campbell
Major Department: Agricultural Engineering

Field experiments were conducted to determine the spatial

variability of the leaching loss of a conservative solute (nitrate or

iodide) in response to a high volume rainfall event. A field plot was

instrumented with tension lysimeters located on 1 m spacing at the nodes

of a 10 m x 10 m grid and two 35 m transects. Irrigation was applied at

4.6 cm/hr for approximately three hours. A chemigation pulse was applied

after the first half hour.

Soil solution samples were collected periodically during

infiltration and redistribution. The velocities of surface-applied solute

pulses (SPV) were determined by fitting an analytical solution to the

convective-dispersive flow equation. Solute mass loss (SML) was estimated

by numerical integration of the soil solution time series and soil water

flux. Soil water flux, determined by irrigation rate, decreased

exponentially with time during redistribution.

Solute mass loss and solute pulse velocities were approximately,

normally distributed. The coefficient of variation was approximately 30

to 60%. The CVs for irrigation and chemigation were 3 to 18%. The

spatial variabilities SPV and SML were determined using regression and

geostatistics. The variability in solute and water application explained

40% of the variability in SML. The semivariogram model explained 20% of

the variability in both variables. All semivariogram models had range

values of approximately 4 m. The nugget values ranged from 40 to 80% of

the variance. There was some indication that reducing the sampling

interval from 1 m to 0.25 m would have reduced the value of the nugget

variance by more than 40%.

Block estimates for SPV and SML using geostatistics are 0-9%

different from estimates using classical statistics. The standard error

of the estimate was reduced 10 to 37% using block kriging.

Uncertainty attributed to intrinsic spatial variability, measured by

the semivariogram, explained 15 to 50% of the variance. Parameter

uncertainty accounted for 5 to 10% of the variance of solute pulse

velocity, and 10 to 32% of the variance of solute mass loss.

A two-dimensional random field with known covariance was simulated

using turning-bands method. The random field was subsampled using

different sample patterns. Experimental semivariograms calculated from

subsamples exhibited up to 40% standard error for small sample size, n -

25. Analysis of cross-validation of fitting the known model to the

subsamples indicated that kriging based on erroneous models was robust to

those errors for fields with smooth semivariograms. Kriging errors

increased by 10% for 1% increases in contamination with outlier values.



Groundwater pollution by chemicals from agriculture on sandy soils is

an ever present problem. Sandy soils are highly permeable and exhibit

poor retentive capability for water and solutes such as nitrate. The

leaching loss of nitrate is important in Florida where frequent, high-

volume rainfall may produce significant leaching and pose a high risk of

groundwater contamination. This is a serious problem where 87% of Florida

households are supplied by public water systems and 94% of dispersed

private rural supplies use groundwater (Heath and Conover, 1981; Leach,


Best management practices (BMPs) have been developed to reduce the

risk of groundwater contamination from agricultural activities such as

chemigation. Best management practices such as drip irrigation reduce

excess irrigation, while plastic mulch reduces evaporation and weed

competition (Bottcher et al., 1986). Plastic mulch reduces competition

for nutrients which may in turn reduce fertilizer requirements. Together,

BMPs reduce the potential for mass loss of applied solutes.

These BMPs often are tested using small plot field experiments

designed to estimate potential leaching losses for various management

practices. Various cultural practices may be compared to identify

practices that produce adequate yields without excessive loss of

chemicals. Selected plots are subjected to intense, high-volume

irrigation representing design storm events likely to result in potential


off-site loss of chemicals. Surface runoff and deep seepage losses are

monitored to estimate actual chemical losses.

Deep seepage losses of chemicals are difficult to monitor. Soil

solution or groundwater samples can be collected to determine temporal and

spatial distribution of the chemicals. In selected locations groundwater

flow can also be measured. Solute transport in the root zone or vadose

zone is difficult to monitor; therefore estimates of solute mass loss must

be numerically simulated. Simulation models range from simple one-

dimensional displacement to complex multidimensional models that account

for hydrodynamic effects, hysteresis, heterogeneous catalysis, and

immobile water for steady-unsteady and saturated-unsaturated flow.

Although the complex models may accurately describe solute transport in

laboratory soil columns, it may not be possible to describe this flow at

the field scale. Simple displacement models, which in the extreme suggest

that solute transport is a simple function of applied irrigation volume,

may be adequate to describe soil water flux for a conservative solute in

a spatially varying media.

The results of the small plot studies may be difficult to interpret

due to the spatial variability of soil properties. High variability

requires a greater number of samples in order to detect significant

differences among treatments. Solute concentrations in the soil profile

are generally variable (Jury et al., 1982) in space and time. In one

study measurements of nitrate concentrations in soil water samples located

at the same depth separated by a few meters were reported to vary by 100%

(Bottcher and Campbell, 1981). It has long been known that hydraulic

properties of the soil are spatially variable. Infiltration rate and

hydraulic conductivity are highly variable (Nielsen et al., 1973).

Spatial variability may also occur in micro-relief of the soil surface

affecting the variability of infiltration. The high degree of variability

of soil properties results in highly uncertain areal estimates of solute

leaching. Where possible, identification and quantification of the

sources of variability should improve estimation procedures.

The variability of the soil properties may not be random. Many

macroscopic properties exhibit properties exhibit spatial structure, or

correlation among neighboring locations (Vieira et al., 1981). The

spatial correlation may result from agricultural practices or internal

soil structure. The range of spatial correlation determines when samples

are independent and can be analyzed using classical statistics (Rao and

Wagenet, 1985). The spatial correlation may provide additional

information to determine the optimal sample size and locations of

observations (Olea, 1985). If significant spatial correlation exists,

geostatistics and kriging can be used to calculate areal estimates

(Journel and Huijbregts, 1978).

The objective of this study is to determine the degree of spatial

variability of solute leaching resulting from high-volume irrigation.

Where spatial structure can be identified, this information would be used

to develop improved estimates of solute leaching losses. These strategies

would provide insight concerning where, when, and how to sample to obtain

optimum information on solute leaching losses.

In this research, the approach was to determine the spatial

variability of solute leaching on a small plot for an irrigation event.

The selection of a simple one-dimensional transport model was to estimate


solute loss that required minimal spatially distributed input data. The

spatial variability and structure of solute mass loss were characterized.

Once the spatial correlation has been modelled, both extrinsic and

intrinsic sources of variability may be identified. Solute mass losses

were evaluated based on spatial structure and correlative relationships

developed with the sources of variability. With this information

alternative sampling strategies were evaluated. Based on the spatial

variability, uncertainty associated with the areal estimates were




Estimation of solute leaching requires an investigation of the

spatial variability of soil properties and the behavior of solute

transport in soil. Soil properties generally tend to be highly variable

and exhibit some degree of spatial correlation. Given the variability of

soil hydraulic and chemical properties it is necessary to select a simple

method for describing solute transport. A simple solute transport model

would depict the spatial variability of solute loss while presenting an

adequate description of water and solute transport. This simplification

is necessary in field sites where it is not practicable to investigate the

spatial variability of soil hydraulic properties.

The spatial dependence of solute leaching may be best described

using regression analysis and geostatistics. Geostatistics provides a

means for optimal areal estimation and for investigation of improved

sampling designs. The uncertainty of the areal estimates is due to model

assumptions, parameter uncertainty, and spatial variability.

Solute Transport

Empirical Results

During an irrigation event, solutes are transported due to soil

water flux. The movement of surface-applied solutes during infiltration

can be estimated by assuming piston-like displacement of initial water by

the invading water (Warrick et al., 1971; Kirda et al., 1973; Ghuman et


al., 1975; Rao et al., 1976). The solute peak will move with the wetting

front in dry soil but, for initially wet soil, lag behind the wetting

front due to displacement of solute free solution initially present in the

soil (Ghuman et al., 1975; Ghuman and Prihar, 1980). Where complete

displacement of initial water occurs, soil columns with either uniform or

depth-dependent water content profiles have identical effluent

concentration distributions (de Smedt and Weirenga, 1978).

The depth and rate of water transport during infiltration depends on

the rate of irrigation, volume of application, and soil water content.

For infiltration, Wood and Davidson (1975) and Bresler et al. (1983)

reported wetting fronts infiltrated to greater depth at lower rate of

water application for equivalent volumes infiltrated. The rate of

application determines the surface water content (Kirda et al., 1973). At

infiltration rates less than the infiltration capacity of the soil, the

surface water content would be equivalent to that associated with a

hydraulic conductivity equal to the wetting rate (Youngs, 1960). Thus

solute will tend to be leached deeper when water is applied to dry soil

than when applied to wet soil (Sommerfeldt and Smith, 1973; Kirda et al.,

1974). The depth to which the solute front is leached depends on the

amount of water added and the average water content behind the wetting

front, but not on initial water content (Warrick et al., 1971; Kirda et

al., 1973; Ghuman et al.,1975; Rao et al., 1976). The antecedent soil

water is displaced or by-passed by the applied water. In either case,

water applied following the solute pulse determined the leaching depth.


For long-term redistribution in sandy soil, irrigation intensity

does not influence the location of the maximum concentration of a

conservative for equal volumes of applied water (Selim et al., 1976).

Similarly, Wierenga (1977) showed that solute movement under transient

water flow is primarily a function of cumulative water input. Kirda et

al. (1974) found this to be true regardless of the initial water content.

Evans and Levin (1969) investigated the effect of different infiltration

rates for the same volume on solute distribution in the soil profile.

They found the final solute distributions after infiltration and

redistribution were equivalent for surface applied solutes when the total

time for irrigation and redistribution was equal. Dahiya et al. (1984)

compared movement of uniformly mixed salt in soil columns for infiltration

and time-matched infiltration-redistribution. Measurement of the center

of mass agreed well with expected movement of a piston front. There was

no significant effect due to water application rate or initial water

content on salt leaching by infiltration followed by redistribution.

Piston type displacement may not occur in all soils. In many soils,

preferential flow through macropores produces incomplete or preferential

displacement of initial water and rapid transport through highly

conductive pore space. This problem may occur in structured or aggregated

soil which exhibits a uniform macroporosity at small scale, or in soils

that have channels due to biotic factors such as roots, worms and rodents.

The macropores in these soils may not be uniformly distributed in a small

area (White, 1985). A fast transport component due to preferential flow

produces a sharp peak at relatively small travel times, whereas the slow

component due to diffusion into immobile water contributes a broad


asymmetric hump extending to large travel times (White et al., 1986). In

sandy soil with less developed structure and minimal aggregation, the two

components produce a single broad peak exhibiting positive skew (Wild and

Babiker, 1976). The depth of penetration of the solute front can be

underestimated without considering preferential flow. Warrick et al.

(1971) reported field measured Cl concentration peaks penetrated deeper

than expected due to flow through large pores.

Convective-Dispersive Transport

Solute transport is determined by the processes of convection,

dispersion and solute interaction. Convective movement of water

transports dissolved solutes through the soil. Dispersion resulting from

hydrodynamic effects and molecular diffusion spreads the boundary of the

solute front. The behavior of solute transport through a soil column can

be described by the process of convective-dispersive flow (CDF) (Warrick

et al., 1971; Ghuman et al., 1975; de Smedt and Wierenga, 1978). The CDF

process has been demonstrated with soil in laboratory columns for steady

(Nielsen and Biggar, 1961) and unsteady (Selim et al., 1976), one-

dimensional horizontal and vertical flow:

R ac = D a2c i ac (1)
at 8z2 az

which describes the change in solute concentration c (mg 1-1), with time

(t) and depth (z) where water content 9 (cm3 cm-3), pore water velocity v

(cm hr-1), and the dispersion coefficient D (cm2 hr-1) remain constant.

This equation assumes that R is a constant independent of C and that local

equilibrium occurs. For conservative solutes, the retardation factor R,

is unity. Convective flow of water is the dominant process. The


dispersion coefficient is a parameter that combines the effects of

processes that cause spreading of the solute front.

Although the CDF equation works well for column studies and has been

used to model experimental transport through soil (Biggar and Nielsen,

1967), there are problems with its use in general field application.

Khanna (1981) reported that the CDF model used for describing movement of

nitrate in homogeneous, nonaggregated soil for steady state soil water

conditions greatly underestimates movement in the field.

Unlike pore water velocity, the dispersion coefficient is difficult

to define and use in a field experiment. The dispersion coefficient is a

term that includes the effects of hydrodynamic dispersion, effective

molecular diffusion and diffusive transport into immobile water

(Passioura, 1971). Hydrodynamic dispersion is the dominant process during

infiltration and the early stage of redistribution. Molecular diffusion

is an important process during low velocity flow and in soil peds.

Laboratory analysis of dispersion indicates that dispersion is

proportional to pore water velocity (Bresler et al., 1983; Davidson et

al., 1983). However, in structured soil or under extreme flow conditions

the dispersion coefficent may not be a simple function of pore water

velocity (Davidson et al., 1983). In steady, saturated flow some

investigators have reported that dispersion is relatively constant,

independent of velocity (Elrick et al., 1979; Smiles et al., 1981). This

may occur in very high flow, near turbulent flow. However, in unsaturated

soil the relationship between dispersion and velocity is not clearly known

(Yule and Gardner, 1978; Smiles et al., 1981; Watson and Jones, 1982).


In field soils the dispersion coefficent has been reported to be

substantially greater than in the laboratory studies, and exhibits both

spatial and temporal variability (Biggar and Nielsen, 1980). In small

field plots, Warrick et al. (1971) found the values of D increased with

time, with solute displacement distance. In the field the dispersion

coefficient varies greatly from location to location (Beese and Weirenga,

1980; Klotz et al., 1980). There is a great increase in dispersion in

field plots, as much as 80-fold over values obtained in laboratory studies

(de Smedt et al., 1986). The increased values of dispersion in the field

are probably due to unconfined lateral flow and the presence of macropores

which increase hydrodynamic dispersion. Dagan and Bresler (1979) found

pore scale dispersion as represented by the dispersion coefficient may be

much smaller than actual solute dispersion encountered on a large field.

The CDF equation applies reasonably well to transport during

transient flow. During transient flow following an irrigation event, the

velocity and dispersion change with time and depth. An effective solute

pulse velocity can be estimated, although values of dispersion are

difficult to interpret. Warrick et al. (1972) described Cl movement in an

undisturbed field using time dependent flow velocity and dispersion and

found good agreement with field data. Warrick et al. (1971) and de Smedt

and Wierenga (1978) showed that the steady state assumption can give a

good approximation to the transient case. Helling and Gish (1986)

determined that the steady state solution also works reasonably well for

cyclic wetting and drying. De Smedt et al. (1986) successfully modeled

transient flow in an experimental sand column. They used an average water

content and water flux in the model. Similarly, Weirenga (1977) modeled


solute movement during transient flow by a steady model using average

water content and showed that peak solute concentrations were greater for

steady state flow, while transient flow produced more spreading.

The CDF equation does not describe well solute transport in soil

which exhibits preferential water flow. The depth of penetration of the

solute front can be underestimated without considering preferential flow.

In structured soil, the CDF model underestimated nitrate leaching in the

early stages of drainage (Barraclough et al., 1983) due to flow in

macropores. The macropores in aggregated soil are generally distributed

in a known manner. Preferential water flow may also occur due to the

presence of root or worm holes that may not be uniformly distributed.

Similarly, in a field lysimeter a simple convective model underpredicted

solute leaching depth (Gish and Jury, 1982).

The CDF equation has been modified to describe solute movement in

structured soils. Water and solute may be partitioned into a mobile

region, where convective-dispersive transport and an immobile phase occur

(van Genuchten and Wierenga, 1976; Gaudet et al., 1977; de Smedt and

Weirenga, 1979). Jaynes et al. (1988) described movement by a simple

convective model assuming a reasonable choice of mobile and immobile water

partitioning. The results of the simulation were sensitive to the

variation of pore water velocity.

An unmodified CDF equation may be adequate for some descriptions of

transport. De Smedt and Weirenga (1978) simulated solute movement using

CDF without accounting for immobile water. Without considering diffusion,

the model still adequately described the location of the solute front.

The study indicated the successful use of convective transport in short-


term transient flow. Under bare soil with essentially no immobile water

the piston-flow model was adequate (Gish and Jury, 1982).

Solute Transport Modelling

Mathematical models are used to simulate solute transport. Since

direct measurement of solute mass loss is impractical, simulation models

are used with field data to calculate solute mass loss. Several

approaches have been developed to provide solutions for various boundary

conditions. Solute transport models based on descriptions of the

processes controlling solute movement have been solved analytically (van

Genuchten and Alves, 1982) and numerically (Bresler and Laufer, 1974;

Burns, 1974; Addiscott, 1977). For unsaturated transient water flow,

these models require solution for water transport in conjunction with

solute transport. These models have been adapted by (Bresler and Dagan,

1983a, 1983b) to provide means for modelling solute transport in spatially

variable media. Stochastic models, which preserve statistical

distributions of the soil properties, have been used to obtain areal

estimates of solute transport. Each of these approaches is discussed


Analytical Solutions

Analytical solutions of the CDF equation have been developed for

several simplified boundary conditions and CDF formulations involving

steady water flow (van Genuchten and Alves, 1982). The analytical

solution for the semi-infinite case with flux boundary conditions is a

commonly used solution. Analytical solutions are applied by fitting the

solution to the time series solute concentration data, calibrating one or

more model parameters using some form of nonlinear least squares (Parker


and van Genuchten, 1983) or maximum likelihood estimation (Knopman and

Voss, 1987). This approach can be used to calculate the time at which the

maximum concentration arrives at a given depth (van Genuchten and Alves,


Analytical solutions have been fit to both resident-averaged and

flux-averaged solute concentration time series. Flux-averaged solute

concentration samples are those collected from actively flowing water.

These samples are typically collected from the discharge of soil columns.

Resident-averaged samples provide a measure of solute concentration in all

water resident in the soil. Resident-averaged samples may be collected by

pore water extraction techniques.

An analytical solution for CDF for flux-averaged concentration (Cf)

presented by van Genuchten and Alves (1982) has been fitted to the

experimental data. This CDF equation is:

Racf = D O2cf v acf (1)
at 8z2 az

subject to the following initial and boundary conditions:

Cf(z,0) = ci (2)

acf(-,t) finite (3)

S o 0 < t < to
cf(O,t) = (4)
0 t > to

The boundary condition (4) represents a flux boundary, a third-type input

boundary condition for solute concentration. The solution to equation 1

subject to 2, 3, and 4 for cf is as follows (van Genuchten and Alves, 1982):

I c+ (co Ci) A(z,t) 0 < t to

Cf(z,t) (5)

ci + (co ci) A(z,t) co A(z,t-to) t > to

A(z,t) erf Rz-vt + t exp(+i'-) erfc Rz-vt (6)
i2(DRt)] I 2(DRt)

This approach works well for saturated or unsaturated steady flow. Bond

(1986) extended this approach by modelling relatively slow, unsteady

unsaturated flow with the CDF assuming 8 and D remain constant using an

analytical solution in a moving coordinate system. Although simulation

using analytical solutions works well for conservative solute and simple

flow conditions, this approach does not handle complex flow fields.

Analytical solutions do not lend themselves to describing nitrogen

behavior and transport during infiltration-redistribution (Wagenet et al.,


Numerical Models

Numerical models have been developed to simulate soil processes too

complicated to be solved analytically. Numerical methods must be employed

for most realistic problems of water flow and solute transport. For

conservative solutes these include CDF models and convection models. Most

commonly, water transport is solved by Richards' equation and solute

transport is modeled by the CDF equations. Bresler and Hanks (1969)

presented a computer model for estimating distribution of chloride in a

soil column using convection. They numerically simulated the simultaneous

transport of water and solute in unsaturated soils with transient

velocities. Kirda et al. (1973) modelled solute transport in a sandy loam


column using Richards' and CDF. The model results fitted the data well

when the dispersion coefficent was accurately estimated. The model was

sensitive to changes in the dispersion coefficent. The model was improved

by the addition of a source term for the solute. Bresler (1973) presented

a model for simultaneous movement of water and conservative solute under

transient unsaturated conditions accounting for diffusion and dispersion.

He included anion exchange in the CDF model which slightly improved the

fit to Warrick et al. (1971) data. Selim and Iskandar (1981) modelled

movement of N under saturated and unsaturated flow. Solute transport and

transformations of nitrogen have been modelled for steady flow with some

success (Misra et al., 1974; Wagenet et al., 1976). Bresler et al. (1979)

included root uptake with the CDF model.

The accuracy of numerical methods depends on choice of model and

boundary conditions, accuracy of numerical techniques, and accuracy and

precision of soil and plant parameters (Bresler et al., 1979). Lack of

agreement with model results occurs due to variation in velocity

distribution during redistribution (de Smedt and Weirenga, 1978) as well

as spatial variability of soil hydraulic properties (Bresler, 1973). The

degree of complexity of numerical simulations is constrained by the

availability of field data for parameter estimation and model

verification. Description of field scale transport is imprecise due to

lack of highly precise or extensive data (Sposito et al., 1986).

As an example, Tillotson and Wagenet (1982) modelled N movement in

field soils following fertilizer application using CDF with source/sink

terms for nitrogen transformations and solved Richards' equation

numerically for water flow. The model used a single valued representation


of the 9-h-K relationship for determining transport. Spatial variability

of the soil hydraulic parameters was not considered. The ability of the

model to predict nitrate movement became progressively worse with time.

Nitrogen transformation coefficients were identified as the probable major

source of error.

Other modelling techniques have been developed to simulate solute

transport when insufficient data are available for soil hydraulic

properties. Simple chromatographic theory has been used to describe water

and conservative tracer movement through soil (Thomas et al., 1978). This

method is less widely used because it requires adjusting the number of

soil layers to fit the dispersion coefficient. Cameron and Wild (1982)

found that the model did not perform well for intense rainfall events, and

that the application of water must be carefully monitored. The numerical

solution was not stable for large differences in water content between

layers. It is used where solute interaction is required and where various

fractions of water mixing are required (Addiscott, 1977). Barraclough et

al. (1983) used this type model to describe nitrate leaching in a layered

soil. This model worked reasonably well for predicting nitrate loss over

the growing season. It appeared that there were macropores that produced

a large drainage volume with little interaction with the soil. These

results are similar to those of Quisenberry and Phillips (1976).

Burns (1975) estimated leaching using a simple model to describe

leaching based on the fraction of applied water leached from the profile

and initial nitrate concentration. The results from his model indicated

that all soil water contributed to leaching and at least over long periods

a simple mass balance could account for solute loss. This model has been


used successfully to describe nitrate leaching from relatively uniform,

unstructured soil without preferred-flow channels.

Piston Displacement

A simplification of the parabolic CDF model is the pure convection

model. This model, used to predict the movement of the solute peak,

neglects the effects of dispersion. The model can be further reduced to

the piston model for description of solute movement over long periods.

The main assumption is the replacement of the actual moisture content

profile by a fictitious uniform one extending from the soil surface to an

equivalent wetting front at depth L(t) with an attempt to satisfy Darcy's

law in the average over the profile and conserve mass. The solution for

the piston model is useful for soils of large variability and for extended

infiltration periods.

Rao et al. (1976) modelled the transport of a conservative

adsorbedd, nondegrading) solute through soil. They did not consider

spreading of the solute front behind the wetting front. Soil water in all

pore sequences was assumed to participate in transport. Thus the solute

front will tend to lag behind the water front if Oi > 0. Initial water was

assumed to completely displace water ahead of the front. During

infiltration the depth of the wetting front (dwf) at a given time t was:

dwf I / (e8-e1), (7)

and the depth of solute front (dsf) was:

dsf = I/ef dwf [ 1 ei/6f], (8)

The solute front depends on the amount of water applied (I) up to time (t)

and the average water content (ef) behind the wetting front at time (t) but

not on the initial water content (el). Similar results were reported by


others (Warrick et al., 1971; Ghuman et al., 1975; Addiscott and Wagenet,

1985a). During redistribution, the depth of the solute front was

determined by the amount of drainable water, (Of-efc) above the solute

front where 8fc is a measure of field capacity. The piston model described

nitrate movement over a long period such as a growing season for a crop.

The piston model was sufficient to describe the movement of the solute

peak under moderate irrigation; however, it was not used to evaluate

nitrate transport for individual events or for intense rainfall events.

Stochastic Models

The extreme variability of hydraulic transport properties (discussed

below) has led to the development of stochastic models (Dagan and Bresler,

1979; Gelhar et al., 1979; Amoozegar-Fard et al., 1982; Bresler and Dagan,

1983a). Two approaches have been developed. One method has been to adopt

a mechanistic model to describe local solute transport allowing for

spatial variability of the model parameters. The second approach is to

ignore the transport mechanisms and simulate the statistical distribution

of solute movement, a transfer function approach.

Dagan and Bresler (1979) assumed that all variations in travel time

in unsaturated soil under uniform water application were due to variations

in 8 only. This approach provided a probabilistic estimate of the solute

distribution providing that the distributions of hydraulic parameters and

application rates are random and independent. In most cases the

distributions were not random and independent, the distributions were

spatially related and thus an integrated solution was not possible.

Bresler and Dagan (1983a, 1983b) successfully adapted the piston

model of soil water flux to a field with spatially varying properties. A


number of assumptions were made in order to use the model. First, it was

assumed the flow is primarily vertical and horizontal flow could be

neglected. Second, soil hydraulic properties were assumed not to vary

along the vertical profile. Third, K(8) relationship was a simple

function, and the randomness of K(9) stemmed from the stochastic nature of

K,. Saturated hydraulic conductivity was assumed to be log-normally

distributed. Water content was initially constant throughout the soil.

The piston flow type profile was assumed at any point in the field during

infiltration and redistribution. The soil column was defined to have

lateral dimensions large enough to include all water flow paths yet

sufficiently small to be negligible relative to the distance separating

two noncorrelated sampling points in the field. Using this water flow

model solute transport was solved using an approximate closed form

solution to the CDF equation. This approach compared well to numerical

solutions for a field of large variability (Bresler and Dagan, 1983b).

Rather than simplifying the water flow equation solutions for solute

transport can be obtained by varying the coefficients of the CDF equation.

The coefficients of the CDF become single realizations of a random

function parameterized by spatial coordinates (Bresler et al., 1983;

Dagan, 1984). The solute concentration at a point in a soil unit, as a

solution of the CDF, then can be viewed as a single realization of a

random field. The quantitative prediction is made from the ensemble

defined over all points.

The CDF model is highly nonlinear and using parameter values based

on field averages produces significantly different estimates than

averaging the results of the CDF model applied to many locations in the


field (Bresler and Dagan, 1979). Persaud et al. (1985) used this method

to conduct a numerical investigation of transport in a heterogeneous field

simulating values of D and v that were log-normally distributed and

correlated. Higher variance of v resulted in large variance of arrival

times. High variance of D was only important at low v. Ensemble averages

were very different from the average of 200 separate curves. Similarly,

Amoozegar-Fard et al. (1982) showed that using average values of pore

water velocity and dispersion coefficient gave poor agreement with field

scale solute movement. Biggar and Nielsen (1976) simulated the movement

of chloride and nitrate using the CDF model. Mean pore water velocity

calculated with the piston model was 20% different from that used in the

CDF model.

Warrick and Nielsen (1980) compared the effect of modelling solute

transport using mean values measured from the Panoche soil and by

conducting a Monte Carlo simulation with values selected for v and D from

the respective pdfs. The mean solute depth from the deterministic case

was shallower. The stochastic simulation produced more spreading of the

solute pulse than simulation of solute transport using mean values for v

and D. Bresler et al. (1979) similarly concluded average values did not

reflect field behavior. Gish and Jury (1982) reported that models using

pore water velocities calculated as the ratio of local water flux to local

water content will underpredict velocity and overpredict travel time.

Other than a Monte Carlo approach there are simpler methods that

incorporate spatial variability of v and D into the CDF model. Rao et al.

(1977) used a Taylor series method in a first-order approximation to

predict the effect of spatial variability on water flux. This method


worked well for short simulations. It used the mean and variance of the

conductivity and a measure of spatial scale to estimate water flux. The

method was quicker and nearly as accurate as the Monte Carlo approach and

could be adopted easily to a new soil. Addiscott and Wagenet (1985b)

described a method for combining measured values and variances to generate

values of K(8) and soil water flux. This was solved by dividing pdfs for

saturated hydraulic conductivity into sections and calculating results

from each section. This method gives results similar to the Monte Carlo

and Taylor series methods.

The second approach to stochastic modelling has been to simulate the

arrival time of solute pulses at some depth in the soil using the transfer

function method (Jury, 1982). Solute flux is given as the convolution of

the input flux distribution and the travel time pdf (TTF) in the soil.

The solute transport in field soil is entirely characterized by

transforming an input function into an output function. The probability

that a surface applied conservative solute will reach depth, L, after a

net amount of water, I, has been applied is as follows:

PL(I) j fL(I') dl' (9)

where fL(I) is the probability density function of the solute peak travel

density function. Since this is a transfer process all of the variations

in chemical transformations and physical processes are assumed to be

unique functions of I. If fL(I) is the average solute concentration at Z=L

in response to a narrow pulse (6) of solute, the input pulse is CIN =

Co6(I). If this is superimposed on the transfer function for an arbitrary


CIN with constant water flux, the solute concentration distribution across

the field at L is as follows:

CL(I) = CIN(I-I') fL(I') dl' (10)

The inlet concentration arrives at Z=L unchanged because dispersion is not

considered. Dagan and Bresler (1979) explicitly included dispersion by

replacing CIN(I-I') with the appropriate solution to the one-dimension CDF

equation with CIN(I) as the upper boundary condition. Jury (1982)

suggested that inclusion of the CDF solution did not quantitatively

improve the solution.

The model can consider spatially variable water application or input

solute concentration flux in terms of a probability density g(i). The

average concentration can be determined for a spatially distributed

solute at any depth, z, can be predicted for the distribution of applied

solute and water:

CL(I) = CIN(I-I') f L/z g(i) ifL(it') di dt' (11)
0 0

The function fL(I) can be calibrated by measuring the distribution of

depths to the peak concentration.

This method is essentially an application of piston displacement.

The model is a representation of the net effect of soil processes, and

modes of solute input on the lifetime of a solute in a unit of natural

soil. If the only loss of solute mass is through the exit surface, then

the lifetime density function can be interpreted as the distribution of

travel times through the soil column.


Jaynes et al. (1988) found that the transfer function can be used if

v and 9 remain constant with each column over time. Jury and Sposito

(1985) compared solute transport using TTF and CDF. The solutions are

different but require a long flow path to separate. The TTF solution is

asymmetric while the CDF solution is Gaussian. The transfer function

approach has been generalized to include physical, chemical and biological

transformation (Jury et al., 1986).

Soil Water Flux

Richards' Equation

Water movement through soil in response to rainfall is an

infiltration process where the water moves downward due to a hydraulic

head gradient. Soil water flux during infiltration and redistribution can

be described by the Richards' equation using gradients of water content

rather than head (Kirkham and Powers, 1972):

ae a D_(19_ aK(8) (12)
at az a 8z az

where: D(9) K/(8a/ah) diffusivity
K(O) hydraulic conductivity

The diffusivity and hydraulic conductivity relationships are

dependent on soil texture and structure. Diffusivity is greater in

structured soil due to the presence of macropores (Hillel, 1980a).

Diffusivity, which is the ratio of the capillary conductivity to the

specific moisture content, depends on the wetting history of the soil.

The hysteresis of D(9) and K(8) seems to parallel the hysteresis of the

moisture characteristic (Kirkham and Powers, 1972). Generally the

distribution of K(8), 8(z,t) and D(8) in space are unknown.


Analytical expressions have been developed for K(8) and D(e) that

depend on the knowledge of residual water content, er, water content at

saturation, 8s, site specific scaling factors and Ks for each location

(Brooks and Corey, 1964; Davidson et al., 1969). The Richards' equation

has been solved for one-dimensional systems assuming all soil properties

except K(O) remain constant and Ks was randomly selected from a log-normal

distribution (Bresler and Dagan, 1983a; Dagan and Bresler, 1983). This

approach is based on the common observation that the hydraulic

conductivity exhibits far greater spatial variability in field soils than

do water content or matric potential (Simmons et al., 1979; Warrick and

Nielsen, 1980; Russo and Bresler, 1982a; Rao et al., 1983).

Another approach is to determine values that effectively represent

average behavior in the field. Dagan and Bresler (1983) derived field-

effective hydraulic properties and showed that, with strong nonlinearity

of K(8) the effective values are only meaningful under restrictive

conditions. In previous work (Bresler and Dagan, 1979; Dagan and Bresler,

1979; Bresler and Dagan, 1981) they assumed the water flow was steady and

later considered unsteady infiltration-redistribution (Dagan and Bresler,

1983). They determined an equivalent homogeneous soil such that the

solution of the flow problem is identical to the average values of a

heterogeneous field. Determination of the effective values required

stochastic models to derive averaged variables.

In either case the use of the Richards' equation required

simplifying assumptions. Either soil hydraulic parameter values were held

constant or effective average values had to be determined. These

techniques require measurements of e(h), D(e), and K(8) for each soil. A


simplified method for determining soil water flux during infiltration and

redistribution could avoid this difficulty.


During early stages of infiltration, soil water flux is dominated by

the diffusivity term. However as infiltration progresses the influence of

gravity becomes more important than D(9) (Kirkham and Powers, 1972). Soil

water flux during infiltration can be simplified by considering the Green-

Ampt assumptions (Kirkham and Powers, 1972). It is assumed that the

thickness of the wetting zone is negligible and throughout the wetted zone

the water content is uniform and constant. Consequently the hydraulic

conductivity in the wetted zone is constant and the flux density is the

same everywhere in the transmission zone. This description has been found

to be effective for sufficiently long infiltration into a sandy soil where

aH/az = -1.


Vertical redistribution is the predominant process by which drainage

occurs. Soil water drains from the wetted soil through a transition zone

into drier soil below. The description of redistribution has been

approximated by two general approaches (Youngs, 1958). In the first case

the water drains uniformly from soil located above the wetting front such

that the draining soil maintains a uniform water content (Shaw, 1927;

Biswas et al., 1966; Staple, 1966; Gardner et al., 1970). This assumes

that evaporation and plant uptake are negligible. The drier soil below is

wetted to the same water content as the draining soil, all soil above the

wetting front is at the same water content. Since at any given time,

water content is uniform with depth, K(O) is constant and the hydraulic

gradient is approximately unity (Black et al., 1969; Davidson et al.,

1969). This case assumes that the soil is nonhysteretic. Dagan and

Bresler (1983) found this to be a reasonable, simplified description of


In the second case of soil water redistribution, hysteresis plays an

important role in soil water flux. Soil water drains through a pronounced

transition zone into drier soil below. The transition zone remains at a

higher water content than the wetted soil due to hysteresis (Bresler et

al., 1983). Water content in the transition zone remains constant until

the soil above drains to apparent field capacity (Youngs, 1958; Biswas et

al., 1966). Soil water flux varies with depth in the profile and is

greatest in the transition zone. The thickness of the transition zone

grows with time, dampening solute variability with depth (Dagan and

Bresler, 1984). In either case of redistribution, water flux is

effectively constant above the wetting front. The change in soil water

flux can be approximated as step changes of constant flux at short

intervals approaching field capacity.

The change in water content and water flux during redistribution has

been described by simple exponential decays. Richards et al. (1956) and

Gardner et al. (1970) reported that water content in the soil profile, 8,

during redistribution after ponding diminishes with time:

S=- a tB. (13)

Where: t time -since beginning of redistribution

a,B empirical parameters

Biswas et al. (1966) showed a decline in water content in the drained zone

and constant water content in the infiltrating, wetting zone. This e


defined a "short-term" field capacity. The soil water flux was constant

throughout the wetted profile over that short period of time. Biswas et

al. (1966) evaluated redistribution in a sandy loam. They found that

water flux during redistribution could be approximated by equation (15).

In q a b In t (14)

where q is soil water flux density and a and b are empirical coefficients.

Similar results have been reported by Leuning and Talsma (1979) in field

soils ranging from loam to clay.

The parameter a is related to the volume of water available for

redistribution and b reflects the rate of soil water loss. Analogous to

exponential decay, the parameter a provides a measure of the water

initially available in the profile for redistribution. This parameter may

be a function of location (x,y) in a spatially varying field and a

function of depth and time. The spatial variability of a is due in part

to varaitions in soil hydraulic properties. It is also a function of

depth, water available for redistribution becomes greater with increasing

depth. It is a function of time in that the applied volume, infiltration

rate, and initial soil water content are time variant.

The exponential decrease of soil water flux with redistribution time

has been verified in an experiment conducted on a profile of Eustis sand

where Davidson et al. (1972) measured water content at 15 cm increments

for several hours following cessation of irrigation. The site was

initially ponded until the profile was almost saturated. Water content

was measured by the neutron scatter method at half-hour increments for

three hours and periodically for two weeks. From these data, the volume

of water percolated from above a selected reference was calculated for


each measurement interval. For the reference depth of 45 cm, the change

in water content is shown in Figure lb and the volumetric flux is shown in

Figure la.

The soil water flux data can be adequately described by an

exponential decrease with time. The relationship underpredicts flux

during the early phase of redistribution, less than 1 hour, and over

predicts flux for later periods. Overall the relationship predicts the

total volume of flux within 5% in the first 12 hours. Near the beginning

of redistribution, within the first ten minutes, soil water flux is

difficult to measure. The discrepancy between the measured water content

and predicted water content could be due to error inherent to measurement

of water content using the neutron probe. At the beginning of

redistribution instability in water flow due to hysteresis reduces the

potential flux. The exponential decrease in soil water flux with time

described redistribution well except immediately following cessation of

irrigation, in which case soil water flux was underpredicted (Sisson et

al., 1980).

It is difficult to model redistribution of soil water due to the

effect of hysteresis. Perreus and Watson (1977) found that accurate

estimation of solute transport requires incorporation of hysteresis.

Soil-water hysteresis tends to result in higher 8 in the upper profile and

the depth of wetting is reduced (Bresler et al., 1983). This retardation

effect is reduced as the wetting and drying hysteresis curves approach

each other. Hysteresis is more pronounced for coarse-textured soils than

for fine-textured soils (Hillel, 1980). Staple (1966) reported that water

content during redistribution in a sandy soil was approximately constant

a) soil water flux





I I t t I I I I I I
1 2 3 4
Duration, t

5 6 7 8

b) Soil water content

45 cm
- -- 60 cm

Figure 1. Soil water during redistribution for Eustis sand following
ponding: a) soil water flux, b) soil water content (Davidson et al., 1972)






Eustis sand
----- q = 2.889 t

-------------------------------- ---------


__ _


below the depth of irrigation penetration. He successfully modelled

redistribution accounting for hysteresis. However, this required careful

measurements of O(h), D(e), and K(8) for each soil.

Sampling Methods

Direct measurement of solute leaching loss is a difficult task.

Transport studies conducted in large tank lysimeters provides a means of

accurately measuring solute mass. In field studies mass loss can be

measured over small areas using or soil water flux meters (Wagenet,

1986). In general indirect methods must be used to estimate water and

solute flux at multiple locations in the field.

Soil Water Flux

Monitoring soil water flux during infiltration and redistribution is

a difficult task. Water flux can be estimated by monitoring changes in

soil water content and water pressure head. Water flux can be calculated

by determining a soil water mass-balance or estimated by using Richards

equation. Solutions of Richards equation require estimation of the soil

moisture characteristic and soil water conductivity-pressure head, K(h).

Estimation of these data requires considerable field monitoring.

Soil water flux through the profile can be determined by calculating

a water mass-balance using water content data. Changes in soil water

content can be accurately measured using single probe neutron

thermalization (Gardner, 1986), single probe neutron attenuation, or dual

probe gamma sorption techniques (Hillel, 1980a, 1980b). This method

integrates a variable volume of soil to determine the water content. The

volume of soil sampled increases as the water content decreases (Bruce and


Luxmore, 1986). Because of the size and expensive of the equipment this

can not be implemented at many sites during a rapidly changing event.

A more common method for determining soil water flux requires

determining the soil-moisture characteristic, e(h), in the field or in the

laboratory on undisturbed cores. In the field, 6(h) is determine by

monitoring water content and pressure head simultaneously. This requires

a plot equipped with tensiometers and assess tubes for a radiation probe

in close proximity (Bruce and Luxmoore, 1986). This method incurrs some

error due to the difference in sampling volumes for water content and

pressure head. The 9(h) may be developed in the laboratory using

undisturbed soil cores (Klute, 1986). This approach is limited by the

inability to obtain soil cores with the structure intact. Further, the

size of the soil core has a strong effect on the results. Soil cores must

be sufficiently large to capture soil structure. Seldom are the results

from laboratory cores equivalent to field data (Hornsby et al., 1983;

Bruce and Luxmoore, 1986). An advantage of the laboratory approach is

that soil cores may be collected following irrigation events at the same

location that pressure head was monitored. This method assumes that the

soil moisture characteristic is nonhysteretic.

In the field, pressure head during unsaturated flow can be monitored

using psychrometers (Rawlins and Campbell, 1986), tensiometers (Cassell

and Klute, 1986), electrical resistance, heat dissipation, or vapor

equilibrium (Campbell and Gee, 1986). Generally, tensiometers provide the

best means of assessing pressure head in the range of solute water flux in

sandy soil. To accurately monitor a spatially variable field cores should

be collected at each sampling location at several depths in order to


provide accurate means of monitoring flux. To convert the measured values

of h to e, tensiometers are installed in the field to monitor changes in

head during the event. The measured values of h must be converted to true

h accounting for 9 and the soil water flux (Elzeftawy and Mansell, 1975).

The result is considerable disruption of the site to collect sufficient

data to quantify the spatial nature of soil water flux.

Estimation of soil water flux requires determination of the

hydraulic conductivity for unsaturated soil. The relationship K(9) can be

resolved in the laboratory (Klute and Dirksen, 1986) or in the field

(Green et al., 1986). In the laboratory, constant head or falling head

methods may be used to determine K(8). Again it is important to obtain

undisturbed soil cores and the quality of results depends on the skill and

knowledge of the experimenter (Klute and Dirksen, 1986). Measurement of

K(9) in the field has the inherent advantage over laboratory methods of

preserving soil structure and incorporating a larger area. The common

approach is to conduct ponded infiltration with drainage experiment, the

instantaneous profile method (Watson, 1966). This requires in situ

measurement of water content and hydraulic head. Due to the inherent

variability of the soil and measurement 9 and h, error in estimation of

K(8) have been estimated in the range 20 to 30% (Fluhler et al., 1976).

Soil Solutes

Evaluation of solute movement requires repeated sampling to

determine the position of the solute pulse and solute concentrations.

These data can be obtained using soil sampling (Alberts et al., 1977) or

soil solution sampling. The former method provides a vertical

distributions of the solute in time while the latter produces a time


series of solute concentration in soil water for the depth of the tension

lysimeter. Soil sampling does provide a means to determine the total

amount of solute at that location, in mobile phase, immobile phase, and

adsorbed. Soil sampling can result in large sample variances due to small

sample sizes (Sisson and Wierenga, 1981; Hawley et al., 1982).

Measurement of water content and solute concentration as obtained from

soil sampling does not provide a good indicator of solute flux,

particularly where K(6) is highly variable (Addiscott and Wagenet, 1985a).

However, soil sampling is adequate for monitoring long-term movement of

solutes (Shaw, 1962; Lund, 1982).

Soil solution sampling is the common method for monitoring solutes

in soil solution. Vacuum extractors, tension lysimeters, and zero-tension

lysimeters are the most common methods for soil solution sampling. The

zero-tension lysimeter, commonly a porous ceramic plate, or trough, is

employed to intercept percolating water and solutes (Dirksen, 1974; Duke

and Haise, 1973; Jordan, 1978). The device is placed in the soil profile

from a trench adjacent to the site. The tension applied to the lysimeter

is adjusted to tension of the surrounding soil matrix (Duke and Haise,

1973). Soil water leaching through the soil profile is captured providing

a direct measure of leaching loss. Zero-tension lysimeters are difficult

to maintain at the correct tension and it remains unclear what fraction of

mobile soil water has been collected. Because soil water moves through

many pore-sequences responding to different tensions simultaneously,

capturing and quantifying soil water and solute transport remains complex.

Tension lysimeters were found to be the only practical means for

collecting soil water samples under field conditions by Alberts et al.

(1977). Soil solution sampling provides a means to collect many samples

at a single location without disturbing the soil profile. Severson and

Grigal (1976) concluded samples are representative of chemistry of soil

water held at the sample tension as applied through the cup. Tension

lysimeters sample flowing water rather than resident water (Parker and Van

Genuchten, 1983). Soil water solute concentrations change little during

sampling which indicates that the solute pulses do not penetrate the small

pores during transport (Warrick et al., 1971). Tension lysimeters sample

unsaturated flow more efficiently than saturated flow. Ceramic cups often

do not accurately sample rapidly leaching water due to preferential flow

(Kissel et al., 1974; Shuford et al., 1977; Shaffer et al., 1979). Solute

samples may differ between samples collected using large cups versus small

cups (Bouma et al., 1982). Large cups may tend to intercept more

macropores and the samples reflect the difference in soil water solute

content between mobile and immobile water. The use of a tension lysimeter

maintained at a fixed suction causes soil water movement unlike that

occurring in the field, regardless of the uniformity of the soil (Cochran

et al., 1970). Rhoades and Oster (1986) and Hansen and Harris (1975)

recommended that standard methods for the use of samplers be adopted which

specify: the use of cups with uniform permeabilities and size, short

uniform sampling intervals, and uniform vacuum design. If the lysimeters

are purged between samples then residual non-adsorbed solutes will not

contaminate the samples (Alberts et al., 1977; Nagpal, 1982).

Spatial Variability of Soil Properties

Soil physical and chemical properties exhibit considerable

variability. The degree of spatial variability depends on the call of


property. Soil properties can be classified into two general types,

capacity or static factors, and intensity factors. Soil capacity

properties such as texture, bulk density, and water content are generally

regarded as normally distributed with low variability (Anderson and

Cassel, 1986). Generally the variation among soil capacity factors is

small, often less than 50%. Bulk density and water content have been

reported to have CV values as low as 20% (Warrick and Amoozegar-Fard,

1979). Mulla (1988) reported CVs for surface temperature greater than 60%

and CV equal to 30 for water content. Reviewing CVs from other

investigations, Warrick and Nielsen (1980) determined average values for

bulk density varied from 7 to 10%, and for porosity they varied from 40 to

45%. Kachanoski et al. (1985) reported CVs from 8% to 30% for surface

horizon properties: density and thickness. Albrecht et al. (1985)

reported CVs ranging from 1% to 14% for water retention.

Intensity factors however such as soil hydraulic properties are

highly variable in field soils. Soil hydraulic properties such as pore

water velocity, hydraulic conductivity, and apparent diffusion coefficient

exhibit highly skewed field-wide frequency distributions best described by

the log-normal distribution with CV values greater than 50% and often

greater than 100% (Nielsen et al., 1973; Biggar and Nielsen, 1976; van de

Pol et al., 1977; Jury et al., 1982). Other investigators also have found

high variability for infiltration rate (Sharma et al., 1980; Vieira et

al., 1981) and saturated hydraulic conductivity (Russo and Bresler, 1980;

Wagenet, 1985). In reviewing several studies, Jury (1985) found that the

CVs for ponded infiltration averaged 194% while the average CV for steady


state unsaturated infiltration was 62%. Albrecht et al. (1985) reported

CVs of 110% to 300% for saturated hydraulic conductivity.

Variability of solute content in soil, capacity factor, also tends

to be large. Quisenberry and Phillips (1976) reported CVs of 50-60% for

solute content in cores collected in the surface 15 cm and CV of 100% for

depths greater than 15 cm from fields on which long term fertilization and

irrigation had occurred. In sandy soil, highly variable nitrate

distributions were reported among profiles for 15 sites (Shaw, 1962).

Selles et al. (1986) reported CV of total nitrogen concentration of 20% in

surface soil. Trangmar et al. (1987) reported CV for exchangeable calcium

equaled 147%. For surface applied conservative solutes, Jury (1985)

reported a range in CV from 60% to 130% for solute concentrations in soil

from the results of 10 studies. Variability in solute concentrations

reflects variability of both soil capacity and intensity factors.

Spatial variability occurs for small plots as well as at the field

scale. In reviewing a number of studies, Beckett and Webster (1971) found

that up to 50% of the variability of some parameters in a field occurred

within a few square meters. The spatial variability affects the

estimation of effective soil properties from measured values.

Two statistical properties of sample populations for soil

properties, the probability density function (pdf) and the coefficient of

variation (CV) indicate the degree to which spatial variability affects

areal estimation. High CV values may also indicate that estimation of

mean values within a specified confidence limit will require a large

number of samples (Jury, 1985). High CV values often indicate that the

pdf is highly skewed and often best described using a log-normal


distribution. Correct identification of the underlying pdf is necessary

to correctly estimate the expected value.

Spatial Structure of Soil Properties

Areal distributions of soil properties often exhibit spatial

structure. Spatial structure can be defined as the description of the

spatial distribution of the soil property (Journel and Huijbregts, 1978).

The spatial structure often can best be summarized by a covariance

function that quantifies the dependence between adjacent sites. The

spatial structure may result from biological and chemical processes within

the soil, climatic variation, or land management activities. The

appropriate structural model may consist of a superposition of two or more

models describing the correlation structure among adjacent values.

It is important to understand and characterize the spatial

variability of soil properties. To accurately conduct statistical tests

the degree of spatial structure must be defined (McBratney and Webster,

1983). Classical statistics assumes independence among samples and does

not account for spatial structure. Further, knowledge of the spatial is

useful in designing sampling strategies (McBratney et al., 1981).

Description of Spatial Structure

A few important characteristics describe spatial structure. First,

the distance beyond which measured values are statistically independent is

the correlation distance, or range. This is also called the integral

scale. A second characteristic, nugget variance, defines the portion of

the variance which is independent of distance. The spatial structure

rarely accounts for the total variance. A additional term, sill,


indicates the portion of variance that is spatially dependent. These

characteristics are illustrated in Figure 2.

Various soil properties exhibit spatial structure. Springer and

Cundy (1987) found a spatial structure with a range of 1.5 m for

infiltration rate and sand content data. Anderson and Cassel (1986)

reported the range of less than 2 m for hydraulic conductivity, soil-water

characteristic, texture, and bulk density. Russo (1984b) reported a

spatial structure range for saturated hydraulic conductivity of

approximately 50 m. He reported that other variables -- crop yield,

electrical conductivity and water content -- also exhibited spatial

structure. Bresler et al. (1984) found sand content exhibited a range of

15-20 m and saturated hydraulic conductivity (Ks) a range of 15-20 m.

The spatial structure may result from large scale external sources

(climate, relief, parent material), management activities, or internal

soil conditions (Bresler and Green, 1982). On a small scale, spatial

structure may be determined by soil aggregation and structure caused by

management activities: irrigation, fertilization, and cultivation

(Albrecht et al., 1985). Lund (1982) indicated that long-term uniform

management tends to decrease the spatial variability. Knighton and James

(1985) reported magnitude of spatial variation and range of soil

phosphorus test values decreased after land leveling. Yost et al. (1982a,

1982b) determined that spatial structure extended 23 to 45 km for several

chemical constituents, appeared to be strongly influenced by climate and

cultural factors. Cameron et al. (1979) found that spatial structure was

related to surface water flux variations caused by micro-relief.

a) components

2.0 -


1.0 -


0 1 2 3 4 5

SI I I I I I I I I I 10 I 11
6 7 8 9 10 11

Lag Distance,

b) semivariogram models

- - Gaussian
------- Exponential
- Linear

Lag Distance, h

Figure 2.

Theoretical semivariogram: a) components, b) semivariogram

variance sill,c

model, \A




Some properties do not exhibit spatial structure, but appear to be

random. Tabor et al. (1985) reported soil nitrate showed no apparent

spatial structure. Such apparent lack of structure may represent either

a random distribution or sampling at the wrong scale.

Soil properties such as water content and tension exhibit time

variant or temporal spatial structure. Saddiq et al. (1985) found that at

high water content was small but the variance may exhibit a high degree of

spatial structure. At intermediate water contents there is a high degree

of variability with low spatial structure (Saddiq et al., 1985). At low

water content there is little variability and little apparent structure.

Greminger et al. (1985) examined the magnitude of the spatial dependence

of the soil-water characteristic and found no structure when the soil was

wet but structure developed as redistribution occurred. Surface

microrelief can be the source of nonuniformity (Solomon, 1984). Stern and

Bresler (1983) found that variability in soil water content decreased

following irrigation for sandy soils but increased for clay soils.

Application Uniformity

Uniformity of water and chemical application to the soil is an

important factor accounting for the variability in solute leaching. Jury

et al. (1982) concluded that high variability of solute leaching was due

to a combination of field properties and characteristics of the water and

solute application methods. Stern and Bresler (1983) evaluated the effect

of nonuniformity of irrigation on soil moisture content and crop yield.

Trends in applied water accounted for approximately 40% of the variation

in yield and water content. Little work has related application

variability to solute mass loss.


Uniformity of application of water and fertilizer is affected by

wind and system design (Solomon, 1979). The system design includes type

of irrigation, spacing of outlets and operating pressure. Wind effect may

be difficult to control and has a significant impact on uniformity.

Vories and von Bernuth (1986) reported that CV increased from 14% to 35%

with increasing wind speed from 1 to 8 m s-1. Bautista and Wallender

(1985) found that uniformity in infiltration also depended on the initial

soil water content, CV of 43% for a relatively dry soil, while CV of 15%

for steady state infiltration.

Small field plot experiments often depend on rainfall simulators for

application. Rainfall simulators generally have high uniformities.

Shelton et al. (1985) constructed a rainfall simulator that applied

irrigation intensities of 7.6 to 16.8 cm h'1 with application CVs of 16%

and 6% respectively. Moore et al. (1983) developed a rainfall simulator

which produced a CV of about 17%. In field practice, an application CV

less than 8% is difficult to achieve.

Spatial Analysis

There are several purposes for evaluating soil spatial variability.

One application is soil mapping or interpolation (Webster, 1985). A

second application is to provide effective parameter values for use in

modelling (Bresler et al. 1983). Spatial analysis is also conducted to

provide an improved method for areal estimation (Burgess and Webster,

1980). Finally spatial analysis may provide a means for identifying and

controlling undesirable sources of variation (Stern and Bresler, 1983;

Russo, 1984b; Morkoc et al., 1985).


Several methods have been used to describe the spatial extent of

soil properties. The simplest methods such as spline interpolation

provide a means to map the soil variability with a given polynomial from

a known set of points. For example, Yost et al. (1982b) used spline

interpolation to map soil nutrient content. Trend surface analysis

provides another means for describing a simple deterministic effect

(Morkoc et al., 1985). The trend surface approach fits a polynomial power

or Fourier series to the entire data set. Agterberg (1974) used a least-

square estimation (LSE) trend surface to derive a simple unbiased

estimator. Walker et al. (1968) used polynomial trends to describe

variation in surface horizon thickness and soil pH related to geomorphic

features. However, minimizing the differences is not suitable for complex

surfaces (Henley, 1981). The method assumes point data, does not account

for finite sample size and is very sensitive to noisy data (Henley, 1981).

The trend surface is very sensitive to the form of sampling pattern

relative to the degree of the polynomial of the trend surface. Trend

analysis does not consider the effect of correlation among data and

consequently is not optimal.

If the sample data exhibit a cyclic behavior the set of observations

can be examined in terms of the spectra in the frequency domain. Byers

and Stephens (1983) used a times series method in frequency domain to

describe the spatial characteristics of hydraulic conductivity and

particle size. Kachanoski et al. (1985) used spectral analysis to

separate stochastic and deterministic parts of signal noise from random

noise. Applied to soil properties, they found that bulk density exhibited


no spatial structure; however, depth and soil mass had significant

autocorrelation and spectral peaks.

Coherency analysis provides a means to assess the linear correlation

between variables as a function of frequency. Greminger et al. (1985)

calculated coherence between the different series in soil-water

characteristic data. This method has an advantage, the correlation

between two variables can be measured on two transects without having to

pair samples or explicitly define phase differences (Kachanoski et al.,

1985). Folorunso and Rolston (1985) used spectral, cospectral, and

coherency analysis to determine cyclic behavior of observations and degree

of linear dependence between soil gas fluxes, soil temperature, and water

content. These techniques require a considerable number of colinear data.

Autocorrelation provides another means for characterizing the

spatial variability of soil properties. This provides a means of

determining the distance at which adjacent samples are correlated.

Vauclin et al. (1982) used a first order autoregressive model to describe

the variation in soil temperature. Vieira et al. (1981) used an

autocorrelation function to describe the spatial distribution in

infiltration rate.


Geostatistics is a set of procedures for determining the spatial

structure and using that information for spatial interpolation or areal

estimation. The intent is to extract from apparent disorder the major

structural characteristics of the regionalized variable (RV) (Matheron,

1963). The procedures consist of structural analysis and optimal, linear

estimation which is called kriging. The structural analysis consists of


calculating an experimental semivariogram, adjusting the results for non-

stationarity, and fitting a model to the experimental semivariogram. The

advantage that kriging provides over other methods of estimation using

spatial analysis is a measure of the error of estimation (Journel and

Huijbregts, 1978).

Geostatistics is based on the theory of regionalized variables. The

random variable of interest possesses certain spatial structure as a

regionalized phenomena. The regionalized variable (RV) is defined by an

observation associated with a sample of specific size, shape, and

orientation known as the sample support. An RV is defined by a

probability density function (pdf) and spatial covariance. The RV is

stationary to the k-th order if the central moments for the joint pdf for

all observations remain equivalent upon translation (Olea, 1975). For

most sample sets it is not possible to determine this beyond second order.

However if the sample pdf can be considered Gaussian and stationary to 2nd

order then the RV is stationary to all orders (Priestley, 1981). Second

order stationarity requires the variance (Var) and covariance (Cov) to be

invariable under translation (Journel and Huijbregts, 1978).

Var(Z(x)) E{ [Z(x)-m(x)]2 ) (15)

Cov(xl,x2) = E( [Z(xl)-m(xl)] [Z(x2-m(x2)] ) (16)

where: xl, x2 -unique values of RV
E( ) -expectation
m(x) -mean value of Z(x).

Often the variance increases with separation distance among observations.

There are seldom sufficient realizations with which to test the degree of

stationarity. This restriction may be overcome by assuming the intrinsic

hypothesis. This hypothesis is that for all vectors defining the distance


between observations, h, the increment [Z(x+h) Z(x)] has a finite

variance that is independent of x. The RV satisfies the intrinsic

hypothesis if its first order differences are stationary in the mean and

in the variance.

E[Z(xi) Z(x2)] m(h) (17)

7(xi,X2) Var{Z(xl) Z(x2)) for all x. (18)

Where the RV is second-order stationary, the semivariogram, y(h), is

related to the covariance, 7(h) Cov(O) Cov(h) (19)

The semivariogram is preferred over the covariance since Cov(x) requires

a known mean and estimation of the mean is often the purpose of the study

(Journel, 1986). The hypothesis of stationarity is an operation decision

validated a posteriori by its efficiency (Journel, 1985).

Additional assumptions concerning ergodicity and homogeneity must be

made in order to use the results from structural analysis for

interpretation of spatial properties. The assumption of ergodicity allows

the statistical ensemble moments of RV to be estimated from volume

averages of a single sample realization. The assumption of statistical

homogeneity, or stationarity in space, means the statistical structure in

space is identical to statistical structure over the ensemble of all space



The first step in structural analysis is to determine an

experimental semivariogram. The semivariogram can be estimated by:

7(h) = 1 n [ z(xi) z(xi+h)]2 (20)
2n(h) i=l

where: n(h) is the number of sample pairs separated by Ih|.


An example of the resulting semivariogram is presented in Figure 2. This

method of calculating the semivariogram is not robust to outliers or

skewed sample distributions. A single value greater than a factor of

three standard deviations from the mean can produce a distinct

semivariogram where none may really exist. Estimation of semivariogram

can be improved by using the Cressie and Hawkins (1980) estimator or

Huber's M-estimator (Dowd, 1984).

There are three components of the semivariogram to be considered in

the model fitting procedure (Figure 2a). The nugget, 7o, is a combination

of measurement and sampling error, and the presence of significant short-

range variation (Journel and Huijbregts, 1978). The range, a, is the

distance beyond which there is no change in 7; i.e., no spatial

correlation. The sill, w, is the sum of the nugget and the between sample

variance which measured by the semivariogram, -y. The value of y7 is the

proportion of variance explained by the model. The sill overestimates the

sample variance due to the spatial covariance among the sampled points.

Model Selection

The second step in structural analysis is selection of the best

model or mathematical description for the semivariogram. There are many

properties and restrictions to the use of semivariogram models. These are

defined and described in detail elsewhere (Olea, 1975; Journel and

Huijbregts, 1978) and are not presented here. The models given in Figure

2b are commonly used.

It is difficult to get a good fit to any model. Several methods

have been defined for selecting the best model and associated parameters.

Least square regression (Vieira et al., 1981; Yost et al., 1982a; Trangmar


et al., 1987), maximum likelihood estimation (MLE) (Kitanidis, 1983), and

eyeball fitting (Saddiq et al., 1985) have been used to match models to

experimental values. Eyeballing the fit appears to be sufficient since

small errors in parameter estimates produce small changes in kriging

errors for spherical and linear models. Fitting the model using nonlinear

LSE adjusted for the number of data pairs fits the model equally weighted

over the range of data (Trangmar et al., 1987). Kriging works well where

the LSE fit of the model to y(h) is good (Myers, et al., 1982), which

indicates that the semivariogram is the main source of spatial structure

and spatial structure is the main source of variance. The best fitting

model should produce the smallest kriging error. The best method for

selecting a semivariogram model is the cross-validation procedure

(Gambolati and Volpi, 1979). Cross-validation can be used on any model to

find the best parameter set. This method is not used as a confirmatory or

hypothesis testing technique such as jackknifing (Davis, 1987). A

semivariogram model is only considered best under the choice of

discrepancy measure, partition set size, predictive function, and number

of models evaluated (Chung, 1984). In this method each observation is

removed from the data set and the point value is estimated from the

remaining n-l values using the kriging procedure. This is repeated for

each sample point providing n-l differences between the kriged and the

known point values, e. The model and parameter selection can be evaluated

by the resulting errors.

Unbiasedness would have the average error, AE:

AE 1/n Z ei 0 (21)

where: n number of observations.


The accuracy of the model in a mean-squared sense should be as follows:

MSE / 1/n Z ei2 minimum (22)

The consistency of the model could be evaluated by the reduced mean-

square error which compares the sample variance with the estimation

variance at each location:


RMSE 1/n e2 1
i-1 -

where: aU2 variance

Kriging errors should be consistent with kriging variances. The kriging

errors should be normally distributed with mean of 0, and variance of 1

and the kriged values should be highly correlated to the observations.

The MSE should be less than the sample variance, in the interval 1 2

J(2/n). The RMSE is very sensitive to deletion of a few points (outliers)

(Myers et al., 1982). If local stationarity does not exist, MSE will be

much greater than 1 and estimation variance much larger (Vieira et al.,

1983). Cross-validation has been used successfully to select the best

semivariogram (Vieira et

1984). Springer and

infiltration data. The

evaluating kriging error

The kriging errors

data. Problems may occui

a small CV. Large nugge

a small CV indicates

semivariograms, where

semivariogram, indicate

al., 1981; Vauclin et al., 1983; Tabor et al.,

Cundy (1987) modelled the semivariogram for

solution was iterative by adjusting range and

s using cross-validation.

depend on the ability of the model to fit the

r in model selection with high nugget variance or

t values will result in large kriging errors and

a small amount of variability. Smooth





smoothly follows the experimental

error (Yost et al., 1982b). Kriging


variance may be more sensitive to changes in the range than in the nugget

(Brooker, 1986). However, the semivariogram behavior near zero is very

important. If the nugget variance is estimated as zero rather than an

appropriate non-zero value, kriging error will be significantly


High cross-validation errors may result from problems with the data.

Data that are not normally distributed will produce biased estimates from

kriging. Genshiemer et al. (1985) found that this was a small problem for

slightly skewed pdfs but depended upon whether the skewness was caused by

a few outliers or a more general behavior. High MSE values may result

from a large minimum lag distance. If there exists appreciable small

scale variance the nugget variance will dominate the MSE term and kriged

estimates will not be much better than other interpolation methods. The

lag distances should be carefully selected based on the scale of the

target phenomena. Poor identification of the range also can cause high

kriging errors. This may occur from insufficiently long lag distances

that confirm the range, or poor selection of neighborhood size for

kriging. Since the model can not be accurately fit to y(h) values where

n(h) < n/2 it is difficult to estimate the range. This is particularly

troublesome with grid data as opposed to transect sampling patterns. When

the data contain more than 10% random noise and fewer than 100 points, it

may not be reasonable to postulate the model properly from the computed

correlation (Chung, 1984).


The third step to the geostatistical process is kriging. Kriging is

a linear, unbiased estimation process for estimating point or block values


(Brooker, 1986). This is an interpolation process for estimating the

value of the RV at an unmeasured location. The best estimate of an RV at

an unmeasured location is the conditional expectation, E:

E[z(xo)(x (xi, ...,z(xn)]. (24)

The expected value of an unmeasured location is the multiple joint

probability function. If the data are multivariate normally distributed

then the kriging estimator is the conditional mean. Often field data are

not multivariate normal and rarely are sufficient data available to assess

the pdf. Linear estimators which only require secondorder moments are

used to estimate values at unsampled locations. As an exact interpolator,

even when less than optimal, the kriging estimator remains the closest

possible linear estimator (Rendu, 1979). The kriged value is the best,

linear, unbiased estimation (BLUE). The estimate of an unsampled location

is a weighted linear combination:

Z*(xo) = Z Ai Z(xi) (25)

that is unbiased if

E[Z*(xo)] E[Z(xo)] (26)

which requires the following condition on the weights,A,

E Ai = 1 (27)

A major advantage of this method is that it provides an estimation error

associated with the kriged values. The variance of the error of

estimation provides a measure of the accuracy of the kriging method. The

linear interpolation method is the optimal estimator if the weights, Ai,

are selected to minimize the estimation variance:

2 Var[z*o zo] 2 AXi(xi,V) y(V,V) -) -Z Aij 7y(xi-xj) (28)
i i J

where: y(xi,V) average semivariogram between points and domain V

-(V,V) average semivariogram of domain V

y(xi-xj) semivariogram among points

This constrained optimization problem the optimal weights, Ai, are obtained

from standard Lagrangian techniques (Lagrangian multiplier, p) by setting

each of n partial derivatives to zero:

a[E([zi z*]2) 2p Z AJ / ai (29)

This procedure provides a system of n+l simultaneous equations with n+l

unknowns which is the "kriging system":

Zi y(xj,xi) + = (xv,xi) i= 1,...,n (30)
and subject to: Z Ai 0 (31)

These equations are solved once for a general sampling pattern and used to

calculate point estimates. With many points these calculations are

extensive. This has been improved by conducting kriging in small

neighborhoods recognizing that all but the closest points are not

considered (Olea, 1985). The resulting values from the optimized solution

are used in the estimation variance to calculate a kriging variance. The

kriging variance for point samples:

aoK = Z Ai 7(Xj,Xi) + p (32)


Block kriging is a linear combination of sample values in or near

the block used for estimating the average block content. While ordinary

kriging is used for point estimation, block kriging is used for areal or

volumetric estimation. The estimation process for block kriging is

essentially the same as that for point kriging (Journel and Huijbregts,

1978; Rendu, 1978).

In the same manner as the kriging equations are developed for

punctual estimation the equations for block kriging are developed from

minimizing the estimation error variance (equation 28). This equation is

also known as the extension variance because it extends the variance of

the point values to the block. The extension variance provides a means to

convert the semivariogram from one block size to another using the

geometric configuration of the data, and block, the structural function of

the semivariogram (Journel and Huijbregts, 1978).

2E 2 Z Ai~(xi, ) 1(v,V) ZE A ij (xi-xj) (33)
i ij

where: y(xi,v) average semivariogram between points and block, u.

-(v,V) average semivariogram between block, v, and domain V.

y(xi-xj) semivariogram among all points.

The kriging variance for block kriging is given in equation (34).

a= Aik iu + f uu(34)
where: -yi average semivariogram between point and block,
7,, average semivariogram for the block (Rendu, 1978).

The estimation variance for block kriging is always less than that of

punctual kriging because the within-block variance is removed from the

error term.


There are several forms of kriging which provide methods for

handling anisotropy, departure from normality, linear drift, and

estimation in the presence of a covariate. The principles for these

methods are fundamentally the same and are described as ordinary kriging.

Other methods such as universal kriging and generalized covariance kriging

allow for estimation in the presence of drift (described below).

Disjunctive kriging and nonparametric kriging can be used where the sample

distribution is significantly non-normal (Rendu, 1980). These alternative

methods provide useful techniques in limited circumstances.


A problem with semivariogram modelling and kriging is the presence

of nonstationarity of the mean (Armstrong, 1984a; 1984b). Nonstationarity

indicates that the expectation of the mean is not constant throughout the

domain. The variation may be defined as a deterministic structural trend

in the data representing some physical process (Volpi and Gambolati,

1978), or some large-scale, low frequency phenomena that cannot be defined

(Starks and Fang, 1982). The phenomena is referred to as drift to avoid

confusion with trends which are well defined surfaces.

Drift affects both nugget and range estimation and may effect the

shape of the semivariogram (Henley, 1981; Starks and Fang, 1982). The

intrinsic hypothesis allows for a drift in the mean which must be removed

in order to determine the correct model. In the presence of nonlinear

drift the shape of the semivariogram is determined in part by phenomena

that have nothing to do with the nature of the spatial structure (Journel,

1986). Gambolati and Volpi (1979) reported that removing the drift

removed most of the anisotropy at short lags. However, in general,


removal of anisotropy is dependent on the process being analyzed. Models

were linear and spherical before drift removal, and linear afterwards. If

drift is only evident at long lags, it will not affect local estimation

(Journel and Huijbregts, 1978; Burgess and Webster, 1980).

Several methods have been used to remove nonstationarity including

filters, regression, or simultaneous estimation of drift and the

semivariogram. One approach to resolving the problem of nonstationarity

is to perform the spatial analysis over some small neighborhood such that

the semivariogram is not influenced by long-range behavior (Journel and

Huijbregts, 1978). The drift may be removed by using a difference filter.

By calculating the semivariogram on differences between locations simple

polynomial drift can be removed without affecting the spatial structure

(Vauclin et al., 1982). Kachanoski et al. (1985) detrended the data by

taking first or second differences. Other filters include median

conditioning and calculation of a moving average (Vauclin et al., 1982).

Vieira et al. (1983) adjusted for drift with a 29-point moving average.

Hamlett et al. (1986) removed drift from soil-water pressure potential

data by subtracting row and column medians from each data point. Cooper

and Istok (1988a) removed drift using sliding neighborhoods.

Regression or trend surface analysis is a second method for removing

drift. Deterministic correlations that have physical basis can be removed

by a spatial trend surface or regression. Generally it is not possible to

generally regress the RV by the coordinates. The fit of a polynomial

expression does not reproduce the measured values of the RV or produce

uncorrelated random error. The estimates of the trend are biased due to

the presence of covariance. Springer and Cundy (1987) detrended


infiltration and sand data using linear regression. The residuals were

analyzed for spatial structure. Yost et al. (1982b) detrended using

polynomial equations. Knighton and James (1985) removed quadratic drift

and fit a linear anisotropic model by nonlinear least squares. Vauclin

et al. (1982) used a first-order autoregression model to describe soil

temperature. They used a fourth-order polynomial to remove the trend.

The semivariogram of residuals had a greater sill and smaller range.

Davidoff et al. (1986) evaluated spatial structure of soil temperature.

They found the trend and subtracted it as a deterministic component. The

resulting residuals were stationary. Detrending significantly decreased

the spatial interdependency.

Drift may be removed by a process that simultaneously determines the

drift function and the semivariogram. One approach, universal kriging, is

used iteratively to identify the drift components and fit semivariogram

models to the drift residuals (Olea, 1975). Only the form of the drift is

applied to the region; coefficients are independently determined for each

neighborhood (Trangmar et al., 1985). The drift is removed from the

neighborhood rather than the entire region. This method is effective

where the drift has a simple structure (Delhomme, 1979). Often drift does

not follow a well defined polynomial structure. Marx et al. (1988) used

universal kriging with linear and quadratic drift for determining the

spatial behavior of soil pH. They found little difference between the two

forms of drift in terms of mean-squared prediction error using cross-


Where the drift function or the semivariogram are complicated, these

functions may be simultaneously identified by utilizing intrinsic random


functions (Delfiner, 1976). Matheron (1973) proposed a method where Z(x)

is viewed as an intrinsic random function (IRF) which could be made

stationary by a process know as incrementing. A kth-order IRF is defined

as a random process that requires a k-order filtering to achieve

stationarity. This method is limited in that only isotropic models are

available. Drift is considered largely indeterminant; linear combinations

of differences effectively filter out the constant. High powers of X

might be filtered out by higher order differences which is precisely the

principle behind the IRF-k (Armstrong, 1984a). The IRF(k) results in

exactly the same system of equations as filtering or universal kriging.

The performance of the automated IRF(k) programs are less than impressive

(Journel, 1986).

Kitanidis (1985) used a minimum variance unbiased quadratic method

to estimate the semivariogram and then estimated drift using weighted

least squares. This method, a simultaneous estimation of drift and

covariance, is also biased and was solved iteratively. In this method an

approximate set of parameter values for the covariance function are

selected and then drift is estimated using weighted least squares. This

method requires assumptions about third and fourth moments that are often

not met for large samples.

Kitanidis (1983) presented polynomial generalized covariance

functions for which parameter selection can be automated. Where there are

significant problems with noisy data or the experimental semivariogram is

highly irregular the model parameter values must be carefully evaluated.

He used maximum likelihood estimation for producing the best estimates of

model parameters. Problems may also occur where the model becomes over-


parameterized to account for drift. The results become unstable.

Kafritsas and Bras (1981) developed an automated iterative approach for

fitting drift polynomials and semivariogram models simultaneously. This

approach eliminated the need to guess at the correct model while

estimating a biased drift function. Kitanidis (1985) found the iterative

approach of Delfiner (1976) and Kafritsas and Bras (1981) was biased and

highly variable.

In comparing methods, Bras and Rodriguez-Iturbe (1985) found neither

detrending or IRF(k) appeared to overestimate or underestimate point

values in point kriging of rainfall from 9 storms of 21 to 29 data points

using either universal kriging or generalized covariance techniques. MSEs

were slightly lower for detrending suggesting more accurate estimation.

The detrending procedure also appeared to give more accurate error

estimations where RMSEs were closer to 1.0. Values of RMSE greater than

1.0 indicate an underestimation of the actual estimation variance. The

detrending procedure ignores errors in the assumed trend.

In addition to dealing with drift, various alternative forms of

kriging have been developed to deal with other data problems. Disjunctive

kriging can be used when the data are not normally distributed. The

empirical pdf is used to develop a transformation to a Gaussian form.

However, disjunctive kriging can not be used if there is significant drift

(Rendu, 1980). Nonparametric kriging is another approach to handling non-

normal data. Where the coefficient of variation is greater than two,

kriging can be done using an indicator function, L(x,z) = (1,0). In this

process L is set to 1 if the value of the RV exceeded a cutoff value, zero

otherwise. For well-behaved phenomena with CV less than one, kriging


works (Journel, 1983). Nonparametric kriging is used for estimating areas

where a given value has been exceeded. Other methods dealing with non-

normal data include using weighted medians or quartile estimators (Henley,


Partitioning Variance

The variability of soil and crop properties can be explained by

various sources. If the sources are independent, the variance may be

partitioned due to individual factors. Morkoc et al. (1985) reported that

three components of the analysis explained yield variation to a line

irrigation source. The deterministic trend was described by a polynomial

and the spatial covariance was described by two components: an

autoregression and a random component. The degree to which a factor

explains the variance of the dependent variable may depend on the area

selected for analysis. Russo (1984b) found that the variability in Ks and

spatial scaling, a, explained 50% of the variability in soil water

potential in an area equal to 50% of the entire field of study.

Accounting for the integral scale in regions of the field with a

relatively small leached zone in the soil profile, 81% of variability in

soil water potential was explained. The variability in soil water

potential and salinity may be responsible for 97% of the variability in

crop growth.

The sources of variability may be defined by trends, correlation

with other properties, or spatial structure. Using multiple regression

hydraulic conductivity and scaling factor were found to explain 25% of the

variability of extractable soil phosphorus and electrical conductivity

(Russo, 1984b). Kachanoski et al. (1985) reported that using a least


squares quadratic surface to describe microrelief accounted for 25% of the

variability in A horizon depth. Yost et al. (1982b) found that polynomial

expressions of trend explained 40% of the variance in soil chemical


Two or more factors may be used to explain the variance on the

dependent variable. Bresler et al. (1984) conducted a statistical

analysis to determine the relative effect on spatial variability of soil

hydraulic conductivity, variation in salt concentration, and soil texture.

He used correlation coefficients to determine principal factors affecting

conductivity. The variation in salinity explained 10% of the variability

in Ks, while sand content variability explained 24% of Ks.


Estimation of solute leaching loss is qualified by the uncertainty

in the estimation process. The sources of uncertainty can be

characterized as two types: intrinsic uncertainty and information

uncertainty (Dettinger and Wilson, 1981). The intrinsic uncertainty

pertains to the natural spatial and temporal variability of the process

being studied. This uncertainty may be quantified but it can not be

reduced. Information uncertainty results from lack of knowledge of the

system. This may result from an incorrect representation of the

fundamental process such as solute transport. It also occurs when the

parameter estimates for the model of the process are not known with

certainty (Cornell, 1972; Burges, 1979). Determining the uncertainty of

areal estimation requires addressing each type.

Analysis of spatial variability provides a means of quantifying the

intrinsic uncertainty. The semivariogram provides a measure of the


intrinsic uncertainty (Dettinger and Wilson, 1981). Uncertainty

associated with estimating deterministic trends in the data also produces

an estimate of intrinsic uncertainty.

The uncertainty of information is a combination of model description

errors and imprecise parameter estimation. Model errors are due to

inaccurate or incomplete description of the system. These errors are

difficult to quantify. They often overlap with intrinsic uncertainty.

Uncertainty due to parameter error can be estimated by determining the

effect of parameter uncertainty on the uncertainty of the results.

Although these sources of uncertainty may not be independent it is often

reasonable to consider these sources independently.

The uncertainty in estimation using a model is determined by the

uncertainty of the parameter values used in the model (Cornell, 1972).

Var[Y] = Var [f(X)] = Var[X] df(X) 2 (35)
I dX I.

The sensitivity of the model to uncertainty in the parameter estimates is
given in equation (36).

ay [ df(X) ax (36)
1 dX m

for parameters Xi, i- 1,...

With a selected model the uncertainty associated with parameter estimation

may be determined by Monte Carlo simulation or first-order uncertainty

analysis. The sensitivities may also be estimated empirically (Tang and

Yen, 1972). Monte Carlo simulation can be conducted by repeated sampling

of the parameter distribution functions and calculation of the resultant

dependent values. This requires probability density functions (pdf) for

each of the parameters (Burges and Lettenmaier, 1975). Defining the pdf

often requires assumptions about the nature of the parameter pdf.


First order analysis provides a second approach to determine the

impact of parameter uncertainty (Benjamin and Cornell, 1970; Cornell,

1972). This analysis provides a means for estimating parameter

uncertainty where only the first- and second-order moments of the random

variables are available. The first order analysis is based on the first-

order terms of a Taylor series expansion or in a system perturbation


Y f(X) f(mx) + (X mx) df + ... (37)
dX mX

This analysis can be used if the relationship is sufficiently linear in

the neighborhood of interest and the coefficient of variance of X is not

large (Benjamin and Cornell, 1970). In the multivariate case:

Y = f(X1, X2, ...) (38)

and the first-order approximation to the variance of Y is:

Var [Y] f M(X) I2 Xi2 (39)

where the parameters, xi, are uncorrelated. The dispersion of the results,

Y, is proportional to the variance of the independent variables (parameter

estimates) and [(af/axi)Im] which is the sensitivity of model f(X) to the

ith parameter xi. The sensitivity of the model results to each parameter

can be determined analytically by differentiating the model with'respect

to each parameter. The sensitivity is calculated by determining the

change in the dependent variable due to a change in the dependent

variable. The sensitivity is absolute when S af/axi, and relative when

S = (af/f) / (8xi/xi). (40)

Sampling Strategies

Optimization of sampling is an important concern in field

experimentation. Optimal design of sampling networks in time (Knopman and

Voss, 1987) and space (Bresler and Green, 1982; Russo, 1984a) may result

in improved accuracy of estimation and in turn reduce the necessary sample

size. Design of the data-collection-network design has the objective of

altering the network configuration until minimum error is achieved (Bras

and Rodriguez-Iturbe, 1985).

Classical methods assume that the value of the RV is determined by

the sum of a regional mean with variation due to location or treatment

effect and an independent random component. The region can be sampled to

minimize the estimation error associated with each treatment source. A

common approach to sample collection is to partition the region based on

expected sources of variance using a hierarchial analysis of variance

(Klusman, 1985). This is a sound approach where the sources are

independent and identifiable.

Where regionalized variables are considered, the model contains a

deterministic component, a stochastic component, and a random component.

The deterministic component may result from identifiable sources or a

regional trend. Understanding trends at the local scale may be important

or may be considered a source of interference to be filtered out. The

stochastic component, considered here to be the semivariogram, requires

sampling to estimate model parameters and form a basis for kriging.

Burgess et al. (1981) recommend reconnaissance surveying to determine the

nature of the semivariogram before a complete survey to determine the

optimum sampling locations.


Sample size may be reduced by using kriging in areal estimation of

soil properties. Where a significant semivariance is encountered in the

data total sample size may be reduced. McBratney and Webster (1983)

compared standard error estimates based on classical statistics and

kriging. The values were 20 to 40% smaller for kriging. This indicated

that half as many samples would be required for the same confidence about

the expected value for the case where the semivariogram accounted for

greater than 75% of the variance. However, it is not possible to predict

the semivariance a priori and the site must be sampled. Based on a prior

sampling improved selection of additional sample sites is possible

(Bresler and Green, 1982; Olea, 1985).

Without a sufficiently large sample size, areal estimation using

geostatistics may not be of advantage. Kriging performed no better than

LSE on a small sample size less than 50 (McBratney and Webster, 1983).

Cooper and Istok (1988b) found that 36 scattered data points were too few

for fitting a semivariogram model to heavy metal and organic carbon data

from a hazardous waste site. Myers et al. (1982) found while analyzing

geochemical data that with 100 sample pairs semivariograms were

undiscernible. However, once density increased to 1000 pairs the

semivariograms were clearly defined.

One major advantage of kriging is the existence of the mean square

error matrix. This provides a means of determining the accuracy of the

data-collection system before making a single observation. Kriging

variance is calculated for every observation in the field (Delhomme,

1978). The variance is reduced by addition of points where the variance

is greatest. This usually occurs along borders or where observations are


least dense (Russo and Bresler, 1982b). If the n(h) data pairs are not

spread evenly about the area, more sophisticated estimates than the equal-

weighted averages should be considered (Omre, 1984). Van Kuilenburg et

al. (1982) found that additional points added to a sparse region markedly

reduces error.

Sampling Designs

Several approaches have been made to improve strategies for sampling

spatially-dependent regionalized variables. Sampling designs have focused

on systematic grids where triangular grids produce the smallest kriging

errors compared to square grids or transects (McBratney et al., 1981;

Olea, 1985). Tabor et al. (1984) found that sampling from a grid provided

better representation than transects. Nugget values were smaller and

resulted in a better defined semivariogram. McBratney and Webster (1983)

obtained 3- to 9-fold gains in sampling efficiency over classical methods

for simple random sampling.

More recently efforts have been made to match the sampling design to

the field behavior of the RV. Trangmar et al. (1985) recommended placing

four closely sampled transects across a coarser grid. A coarse grid

provides a basis for kriging interpolation and the transects enable

detection of anisotropy and establishment of the semivariogram at short

lags. Where geometric anisotropy occurs an elongated grid should be used.

For block kriging, the optimum sampling pattern is related to the

density of the observations and the block size (Burgess and Webster,

1980). There should be a large number of evenly distributed sample points

within the block for estimation of the block variance. There should also

be sufficient external sample points for estimation of the domain


variance. Where the samples are sparse the estimation variance decreases

with block size until a large number of samples are within the block as

opposed to those surrounding the block. The maximum estimation variance

occurs at the center of the grid cells.

Bresler and Green (1982) developed a method that iteratively adds

sampling points to an existing sampling pattern that maximizes the number

of distance classes and attempts to obtain uniform class sizes. The focus

on class size is drawn from investigations of spectral analysis where

definition of the spectral density function requires careful estimation at

all sampling frequencies. The process begins with a user-determined

pattern based on an attempt to uniformly monitor the field of interest.

The results of their method improved the kriging errors about 20%. This

approach results in a sampling network that is very close to a systematic

network (Russo and Bresler, 1982b).

The iterative selection approach was extended to provide a method

for designing the initial sampling pattern by Russo (1984a). Russo

developed a method for determining the sample locations based on a

constant number of sample pairs within each lag class and uniformity of h

within each class. The method produced smoother semivariograms than a

standard square grid. Warrick and Myers (1987) extended the method to

provide a pattern that was optimized with respect to the distribution of

sample couples within a class size. This provides a means of determining

the class sizes of interest and balancing the sampling to have even


Hughes and Lettenmaier (1981) used an iterative algorithm to adjust

the location of sampling sites in order to minimize regional kriging


variances. The efficiency of their approach depends on the assumed

initial locations. The final pattern depends on the shape of the

semivariogram model. They tested their method using a simulated data set

which possessed a known two-dimensional covariance. Veneziano and

Kitanidis (1982) presented a method for sequential sample location

selection based on minimizing the estimation variance.

Other approaches have been developed that optimize the sampling

scheme according to additional criteria. Using a multicriteria decision-

making function, Bardossy and Bogardi (1983) presented a method for

selecting sampling locations based on spatial dependence and sampling

cost. Where the appropriate model is known, a sampling pattern can be

developed to minimize the estimation variance and the number of observa-

tions. This analysis requires prior knowledge of the semivariogram. This

procedure provided a means to include economic cost and observational

capacity in the site selection process. This was particularly useful for

large sampling programs.

Although the prior discussion concentrates on areal sampling

schemes, temporal strategies are important for estimating solute leaching.

Knopman and Voss (1987) conducted a sensitivity analysis which indicated

the necessity of collecting soil solution samples near the arrival of the

peak to provide the best estimate of pore water velocity. As previously

discussed solute leaching is spatially variable with a lognormal pdf for

peak arrival time. Field scale sampling designs need to consider

irregular sampling schedules.


A common technique for evaluating alternative sampling patterns is

to perform multiple sampling on an artificial data set. This allows for

greater variety of sampling patterns with a large number of points. The

statistical characteristics of the simulated data can be clearly defined

a priori. In particular, a two-dimensional field on known covariance can

be generated. This can be accomplished using the turning-bands method

(Journel and Huijbregts, 1978; Mantoglou and Wilson, 1982a; 1982b). With

this method a unidirectional line covariance is created using a simple

spectral method. The line process is generated discretely. A finite

number of lines are generated radially about a selected reference point.

The two-dimensional field is then created by summing orthogonal

projections from each line to a point on the grid to be created. This is

repeated for every point. The simulation error associated with the

generation of the two-dimensional field due to a finite number of lines

and discretization in the spectral process is less than 2.6% of the sample

variance (Mantoglou and Wilson, 1982b). The result is a field of points

with known mean, variance, and covariance. By sampling the simulated

field and testing the estimation of the statistical parameters different

patterns and methods can be evaluated (Hughes and Lettenmaier, 1981).


It is clear from the literature that soil properties exhibit spatial

structure. However, there have been few studies describing the spatial

structure of solute leaching in sandy soils. Similarly there has been

little work evaluating the source of spatial variability; intrinsic and

extrinsic factors. In the evaluation of alternative management practices,


extrinsic sources, particularly variability in application rate may be


Solute transport during and immediately following irrigation events

is primarily controlled by the convective flow of water. Experience

indicates the conservative solutes travel with the wetting front during

infiltration and redistribution in homogenous soil. Solute transport is

primarily convective and the effective pore water velocity can be

adequately described by the convective-dispersive flow model for transient


Soil water flux during infiltration and redistribution is difficult

to quantify a many locations simultaneously. Although water flux can not

be measured directly accurate, indirect methods provide measurements of e,

K(e), and h(E) that can be used to calculate flux. For an accurate

analysis of spatial variability these relationships must be developed for

many locations in conjunction with measures of solute transport. This is

extremely time consuming, destructive, and requires many soil samples.

This approach is not appropriate for evaluating solute leaching at a large

number of locations.

Soil water flux can be estimated using simplifying assumptions.

During infiltration that water flux is equivalent to the irrigation rate.

And it can be assumed that soil water flux follows an exponential'

reduction with time. In this way an approximate model for soil water flux

can be used with a time series of solute concentrations in soil solution

to estimate solute flux. Solute transport can be monitored at many

locations using soil solution samples for evaluation of the spatial

variability of solute leaching.


The spatial structure can be described several ways. However,

geostatistics provides unique advantages; a measure of uncertainty and

optimal behavior for estimation. The primary tool of geostatistics, the

semivariogram is used to describe spatial covariance. The semivariogram

is likely to be confounded by drift, the presence of a non-constant mean.

There are several techniques available for handling drift; differencing to

adjust for trends, removing trend surfaces from the data, and correlation

against soil or site characteristics. The latter approach has the

potential advantage of identifying meaningful physical relationships.

The semivariogram is used for areal estimation through kriging. A

key element to kriging is fitting a model to the experimental

semivariogram. Although several methods have been developed to fit

models, with or without drift, it appears that cross-validation is the

best available technique. Rather than optimizing the model fit to all

data, cross-validation tests the robustness of the model by removing data

points from the analysis.

Several sampling strategies have been developed to assess spatial

variability. These strategies optimize different factors that influence

the development of semivariograms; number of sample pairs, inter-sample

distance classes, and coverage of the field. The simplest strategy for a

small sample size is to consider a simple grid or transect. Different

intersample distances are also important to consider.



This project was conducted to evaluate the spatial variability of

nitrate and iodide leaching losses in a small field plot during single

irrigation events and to provide areal estimates of solute leaching loss.

The approach in this study was to describe the solute pulse movement

through time at a single depth at many locations. It was assumed that the

solute concentration time series resulted from an integration of the

spatial variability of soil hydraulic and chemical properties for each

location. Two metrics, solute pulse velocity and an estimation of the

mass of solute loss, were selected for evaluation of the spatial

variability solute leaching.

Three experiments were developed to evaluate the spatial variability

of solute leaching loss. These experiments were conducted to determine

the impact of different sampling strategies for similar leaching events.

The first experiment was designed to determine the effect of intersample

distance on the estimation spatial covariance. Three intersample

distances, 0.25, 0.5, and 1.0 m, were used to determine the effect of

intersample distance. In the second experiment spatial variability of

solute leaching was ascertained using a grid sampling format. In the

third experiment spatial variability was ascertained using long transects.

This provided a greater number of sample pairs than grid sampling for

improved estimation of the spatial covariance. Iodide and nitrate were


used as tracers because they behave as conservative solutes and should

reflect the transport of highly mobile solutes.

The estimation of solute loss was based primarily on soil water

samples collected using tension lysimeters during each event. Two metrics

were used to evaluate the spatial variability of solute loss; solute pulse

velocity (SPV), and solute mass loss (SML). The solute pulse velocity was

approximated by an estimate of the effective pore water velocity. An

analytical solution of the convective-dispersive flow (CDF) equation was

fit to the time series of solute concentration from each location using a

nonlinear least-squares fitting technique (Parker and Van Genuchten, 1983)

to estimate effective pore water velocity.

The estimates of the mass of solute leached at each location were

estimated by integrating solute concentration data from each location with

a simple soil water flux model. The soil water flux model assumed flux

was equal to the irrigation rate during infiltration, and estimated by an

exponential decay with time during redistribution. Estimates of solute

mass flux were estimated for each sampling time during each irrigation

event. At each sample location the instantaneous mass flux estimates were

integrated over the time of sampling to determine the mass of solute loss

during the event.

Geostatistical procedures were used to identify spatial structure

and calculate areal estimates of SPV and SML. Simple linear regression

with application variables was used to remove drift from SPV and SML.

Geostatistical methods were also used to conduct statistical simulations

to evaluate sampling patterns and densities used for collecting mass loss

data, and to provide insight in the selection of more efficient sampling.


Each of the procedures used in the above analysis is described below.

The first section describes the evaluation of the chemical analysis for

iodide and nitrate, and iodide and nitrate adsorption experiments. The

next section describes field experiments for evaluating tension lysimeter

behavior and spatial variability of solute movement. Next, the procedures

for estimating solute mass loss during infiltration and redistribution is

presented. Finally, the geostatistical methods are presented.


Nitrate and Iodide Analysis

Soil solution samples were analyzed for nitrate and iodide using ion

chromatography. A Dionex ion chromatograph (Model QIC) with an anion

separator column, fiber suppressor column and electrical conductivity cell

detector was used to determine nitrate and iodide concentrations.

Complete description of this procedure is found in Appendix A.

Solute Adsorption on Soil

In these experiments nitrate and iodide were treated as nonreactive

solutes. To verify this assumed behavior, soil adsorption of iodide and

nitrate from solutions containing 5, 10, and 20 mg iodide or nitrate 1~1

tap water was measured. This experiment was conducted using the double-

centrifuge technique (Dao and Lavy, 1978). Adsorption was determined for

4- and 24- hour periods with 1.5 ml solution in 5.0 g fresh top soil.

From this experiment adsorption was calculated assuming linear adsorption

for iodide and nitrate. The resulting retardation factors were found to

be 1.05 and 1.08 for iodide and nitrate, respectively, thus it was

concluded both ions were essentially nonadsorbed.

Field Experimentation

Several experiments were conducted on a small field plot on the

University of Florida campus. The first set of experiments was designed

to evaluate the use of tension lysimeters, evacuated porous ceramic-cup

samplers, for collecting soil water samples during irrigation events and

determine the best procedure for using them. These experiments and

results are described in Appendix B. The second set of experiments was

designed to determine the spatial variability of solute leaching.

Field Plot Development

The field plot used in these experiments consisted of a 50 by 20 m

bahia grass field located at the northeast corner of Hull Rd. and Museum

Rd. on the campus of the University of Florida. The site was selected as

a one typical of agricultural activity of shallow sandy soil in North

Florida. The grass cover was estimated by University field crew to be

approximately 10 years old. The microrelief consisted of scattered grass

clumps and bare sand. The clumps had an average spacing of 12 cm and a

height of 2 cm. The mosaic of grass and sand was uniform across the site.

The rooting depth of the grass was approximately 50 cm. At an early time

the site had been planted to peanuts. Prior to this study the site had

been occasionally used as a tractor practice driving range by the

Agricultural Engineering Department. No attempt was made to mitigate any

residual effects of prior agricultural activities.

This site was relatively undisturbed. The residual effects of

cultivation and chemigation were expected to be minor. Preliminary

assessment of infiltration rate, physical characteristics, and chemical

properties did not produce evidence of spatial trends in soil properties


on the site. This site should show the response of a reasonably

undisturbed agricultural field and consequently the spatial variability

due to intrinsic sources should be great as possible. The effects of any

cultivation were expected to be a dominant factor and the analysis of

solute leaching loss could be dominated by cultivation. The effects of

specific practices should be investigated in a separate experimental


The underlying soil is a Kendrick sand (Loamy, siliceous,

hyperthermic, Arenic, Paleudult). Several soil samples were collected at

the site to describe selected physical characteristics. Immediately

before and after each solute leaching experiment soil core samples were

collected for moisture content analysis using a thin-walled aluminum tube

and mallet. The average value (Ave) and standard error (SE) were reported

for each property:

n n
Ave = Z Xi SE (Z Xi) / n i = ..., n (41)
i i

The sample holes were refilled with sand. Three holes were excavated at

selected locations after the leaching experiments were completed. From

these holes, four undisturbed samples were collected from each of five,

10-cm depth increments. From these samples the following analyses were

conducted: bulk density by core method (Blake and Hartege, 1986) (Table

1); particle-size by hydrometer and sieve method (Gee and Bauder, 1986);

and soil moisture characteristic (Klute, 1965) (Table 2). This soil

exhibited no significant change in texture or bulk density in the upper

0.8 m indicating the lag of a confining layer. The soil was considered

vertically homogeneous. A limited number of samples were collected for

physical analysis due to the destructive nature of sampling.

Table 1. Average value (Ave) and standard error (SE) for initial water
content (9O), and residual water content (e8), bulk density for field plot.

Initial Water Residual Water
----- cm3 cm3 -
Ave. SE Ave. SE

Bulk Density
- g cm-3
Ave. SE



15.1 0.32
14.2 0.25



Table 2. Soil physical properties of typical Kendrick sand (Loamy,
siliceous, hyperthermic, Arenic, Paleudult).

Particle Size Distribution (% < 2mm)


(dia. 2.0-
mm) 1.0





0 50 (%) 1.6 2.8 25.7 61.2 0.68
(SE) 0.11 0.21 2.5 1.2 0.04

Sand Silt Clay Texture
2.0- 0.05- <
0.05 0.002 0.002

93.7 0.35 5.95 fine
1.1 0.03 0.50 sand

Soil moisture characteristic curve

head (cm)

0 20 (cm3 cm-3)

3.5 20 30 45

38.0 33.6 30.5 2
1.2 1.1 0.6

60 80 150 200

5.1 20.3 17.1 13.3
0.7 0.5 0.25 0.15

- cm -


Field Leaching Experiments

Three experiments were conducted during 1985-1986 to evaluate the

spatial variability of nitrate or iodide leaching loss during and

immediately following high intensity irrigation on a sandy soil. In each

experiment the placement of tension lysimeters was different. Experiment

I was conducted to evaluate the effect of small inter-sample distance on

sample collection and spatial structure. Experiments II and III were

conducted to evaluate the spatial variability of solute leaching and the

use of geostatistics for areal estimation. Data from the short transects

were used to evaluate the spatial structure at shorter lag distances.

The three experiments consisted of different placement patterns of

lysimeters sequentially in a single field (Figure 3a). Experiment I was

conducted on a single 5 m transect with lysimeters placed at 0.25 m

spacing (Figure 3b). In Experiment II the lysimeters were placed in a

square grid 10 m x 10 m with 1.0 m spacing between nodes (figure 3c). In

Experiment III the lysimeters were placed in two parallel 35 m transects

(separated by 2 m) at 1 m spacing (Figure 3d). Ten additional lysimeters

were placed in a portion of the north transect to produce a 12.5 m

transect with intersample distance of 0.5 m.

The experimental conditions for each experiment were similar (Table

3). The experiment was designed to evaluate the potential nutrient

leaching of a conservative solute applied as chemigation and followed

immediately by a high-volume rain. This would produce solute leaching

under one extreme scenario. In these experiments an initial 45 minute

irrigation was applied to wet the soil before each experiment.

Immediately following the prewetting, solutions of sodium nitrate or


potassium iodide were applied in a short pulse of chemigation through the

irrigation system. Following chemigation the site was irrigated with

several centimeters. Soil solution was sampled during irrigation and


Irrigation and Chemigation

The irrigation system used in this study was set up to provide

approximately the volume of water that would fall during a 2-hour, 10-year

return period storm (Weiss, 1962). The method of Weiss (1962) as

described by Schwab et al. (1981), was used to calculate the necessary

irrigation depth of 9.7 cm for such a storm. Irrigation uniformity was

controlled by location and number of sprinklers, and applied pressure.

Water was applied through impact sprinklers (Rainbird model 25ASFP-TNT

with 0.32 cm nozzle) located along the sides of the plots outside a 2-m

buffer strip. The sprinklers were clustered in groups of four on top of

five-foot risers. The water supply was received through a 5-cm branch

from the City of Gainesville water main at approximately 50 psi (although

the city main pressure was greater than 80 psi, head at the site rarely

exceeded 50 psi). The line pressure varied somewhat, approximately 5 psi,

during each experiment and adjustment of a valve compensated for pressure

variation. The location of the sprinklers for each experiment was

determined by the size and shape of areal distribution of lysimeters for

each experiment, and the application distributions of the available

sprinkler heads. The greatest uniformity was achieved with four heads per

side located at 16-m spacing with lateral pressure of 40 psi for each

experiment. Uniformity was determined for each leaching experiment by

a) Areal view of lysimeter placement for all experiments.

o 0

0 0 00 0 a 0

1 8


o o o
0 0 O


** ** *

For Detail See Inset

* o Lysimeters, 1.0 m Spacing

b) Areal view of lysimeter placement for Experiment I.


0 meters

0 0 0 0 0

0ooooooooooo o olololololololololo o 0 o

0 0 0 0 0

8 20


o Vertical Lysimeter, 50 cm Depth
| Lysimeters Placed at 45 From Vertical, 50 cm Depth
(see inset)

o 25 Lysimeters, 0.5 m Spacing

olo 20 Lysimeters, 0.25 m Spacing

Figure 3. Experimental site layout. a) Areal view of lysimeter placement
for all experiments. b) Areal view of lysimeter placement for Experiment
I. c) Areal of lysimeter placement for Experiment II, d) Areal of
lysimeter placement for Experiment III.

o 0 0
o 0 o
0 0 0
o 0 0
0 o o
0 0 0
0 0 0


c) Areal of lysimeter placement for Experiment II

Experiment II, Square Grid


o Vertical Lysimeters at 1.0 m Spacing

d) Areal of lysimeter placement for Experiment III.

Experiment III, Long Transect



0000000 000



Figure 3. Continued

09 0 0 0
.0 0 0 0 0





0 0 0

0 0



0 0 0




calculating the coefficient of variation for application volume collected

near each of the tension lysimeters.

Nitrate and iodide solutions were applied by chemigation through the

irrigation line. The solute solutions were mixed in concentrated

solutions and injected into the main irrigation pipe sufficiently upstream

to insure adequate mixing. Actual application rates were determined by

collecting irrigation water samples from the plot. The actual mass of

solute and irrigation volume delivered to the site during each experiment

was determined from 250-ml open-top collection bottles placed near each

tension lysimeter prior to irrigation. This method was used in all

experiments except Experiment B8 (Appendix B).

In B8 a combined solution of iodide and nitrate was applied to the

plot using a hand sprayer at an areal application rate of 40 kg ha-1.

Filter paper randomly placed on the small (2.5 m x 5 m) plot were analyzed

to determine the uniformity of application. Although the hand sprayer was

an effective means for uniform application on a test plot it was grossly

inefficient for larger areas.

A summary of the conditions of each experiment are presented in

Table 3. A detailed discussion of the variability of application is

presented with the discussion of the spatial variability of solute

transport in the next chapter.

Soil Solution Sampling

Tension lysimeters were used to collect soil solution samples

simultaneously at each node during each experiment. Vacuum was applied

through a manifold to each lysimeter. The lysimeters were fabricated by

cementing porous ceramic cups with dimensions of 5-cm long and 0.2-cm o.d.

Table 3. Irrigation and chemigation parameters for field experiments.

Exp. 1 Exp. 2 Exp. 3

Number of 20 98 96
sampling sites

Irrigation hr 3.83 4.58 3.83

Chemigation hr 0.62 0.33 0.33

Post-Chemigation hr 2.34 3.91 2.83

Average cm hr-1 4.60 4.61 3.58
Irrigation Rate

Solute Applied Nitrate Iodide Nitrate

Average Areal kg ha-2 5.61 12.2 5.15
Chemigation Rate

Solute Concentration mg 1-1 17.7 80.3 43.6
in Chemigation

Number of samples 17 12 20
Collected per site

with nominal 1-bar bubbling pressure (Soil Moisture Equipment, Inc.) to an

open end of a 70-cm length of 0.20-cm o.d. plexiglass tubing. The

diameter of the cups was approximately the same as the tubing outer

diameter. The upper end of the tubing was sealed with a plexiglass plate

through which pressure fittings were placed. Nylon tubing (0.05-cm o.d.)

was placed through these pressure fittings extending down the plexiglass

tube to the bottom of the ceramic cups. The other end of each nylon tube

was connected to a sampling manifold located 2 to 6 m from the lysimeter,

depending on the sample location within the plot (Figure 4).


Each sampling manifold consisted of a 50-cm section of 3 cm i.d. PVC

pipe sealed at both ends. The pipe was attached to an electrical vacuum

pump using 0.64 cm i.d. nylon tube through a pressure fitting. Brass

nipples were attached along each pipe and a short length of tygon tubing

was attached to the nipple. A pinch clamp was placed on each piece of

tygon tube. A short piece of glass tube with a neoprene stopper was

attached to the tygon tube. The nylon tubing from a lysimeter was placed

through the stopper and extended 3 cm beyond the end of the glass tube.

A 50-ml polycarbonate centrifuge tube was firmly attached to the stopper.

The polycarbonate centrifuge tube was selected for volume, shape of

orifice (ensuring a secure fit with neoprene stopper), and non-reacting

nature of material. Vacuum was applied to the lysimeter through the

centrifuge tube. The centrifuge tube was removed momentarily to release

suction. The pinch clamp was used to isolate each lysimeter from the

vacuum system. Two or Three sampling manifolds were used to form a

station. Four stations were established outside a 2-m buffer along the

sides of each experimental plot for a total of 10 manifolds servicing

approximately 100 lysimeters.

A tension lysimeter was installed at each node as indicated in the

sampling pattern for each experiment and was located at a depth of 0.5 m

(Figure 4). In Experiment I (Figure 3b) half of the lysimeters were

inserted at a 45" angle with the ceramic cup located at 50 cm depth so as

to test the impact of vertical placement (see Appendix B). In all other

experiments lysimeters were inserted vertically.

Care was taken during installation to minimize soil disturbance. At

each node a hole of slightly smaller diameter than a lysimeter was punched

To Pump

1/4 Nylon

3 5m x .05cm Nylon Tube

. (9) Additional
Collection Tubes


Sampling manifold and collection tubes for tension lysimeters.




o:n, 11 4

:' : .: '


Figure 4.


using a thin-walled aluminum tube and mallet. The lysimeter was pushed to

the bottom of the hole resulting in a tight fit. Loose soil from the hole

was formed in a small mound around the neck of the lysimeter to prevent

leakage down the tube. Although maintenance of the site required foot

traffic across the site there was no evidence during irrigation of

compaction or soil sealing problems. Near the sampling manifolds where

considerable traffic occurred there was ponding during irrigation events.

During each event, samples were collected at irregular intervals

during the last part of irrigation and during the major phase of

redistribution (Table 4). After approximately ten hours from cessation of

irrigation the soil water content became too low to collect representative


Table 4. Sample collection times for each experiment in hours past
cessation of chemigation.

- - - - - - hours - - - - - - - -

Exp I 0.50 1.00 1.40 1.90 2.50 3.00 3.60 4.25 4.75 5.33 5.75
6.25 6.75 7.25 9.25 11.00

Exp II .67 1.58 2.42 3.08 3.67 4.33 5.0 5.67 6.33 7.08 7.83 8.83

Exp III 2.5 3.00 3.33 3.77 4.00 4.33 5.00 5.67 6.33 7.08 7.50
8.0 8.58 9.33 10.08 10.83

At each sample collection time, vacuum was applied to the vacuum

system simultaneously at all lysimeters for approximately one minute.

Suction was released and immediately re-applied. This resulted in an

initial purge of water from the cups which was discarded. Vacuum was re-


applied for approximately 10-15 minutes. Every five minutes each sampler

was checked by releasing and immediately reapplying vacuum. This resulted

in collection of between 25 and 40 ml of soil solution. The time of

sampling was recorded as the midpoint of the sampling period. Once

sufficient sample was collected the lysimeter was disconnected from


Each sample was collected into a polycarbonate centrifuge tube.

Samples were transferred to scintillation vials and placed on ice.

Samples were taken to the laboratory and prepared for analysis immediately

following collection.


Tensiometers were used during Experiment II to track the arrival of

the wetting front and the change in tension during redistribution.

Tensiometers were constructed in the same manner as the tension lysimeters

except the 0.05-cm o.d. nylon tubing terminated in a mercury manometer

rather than the sample collection manifold. Unfortunately many

tensiometers did not function well: there was a significant lag time

between arrival of the wetting front and apparent change in soil matrix

potential. There was a high degree of variability in time between when

the tension lysimeters were able to collect the first samples and when the

tensiometers indicated arrival of the wetting front. More importantly,

the lag-time problem was not consistent among tensiometers. Given the

high degree of variability among the tensiometers the data were considered

suspect for these short experiments and not reported.

Solute Leaching Loss

Solute Pulse Velocity

An important measure of the spatial variability of solute leaching

is the variability of the displacement of the solute pulse through the

soil. Weirenga (1977) found that the displacement of the solute pulse

during transient, unsaturated flow could be described by the convective-

dispersive transport equation (equation 1). The estimate of pore water

velocity resulting from that description represents an effective velocity

averaged over infiltration and redistribution. The effective pore water

velocity was used here to define the velocity of the solute pulse.

Solute pulse velocity was calculated for each location for each

experiments. The SPV values were determined at each location by fitting

an analytical solution (van Genuchten and Alves, 1982) of the CDF equation

for flux-averaged solute concentrations to the time series data, c(t),

(equations 5 to 6). This approach provided a means by which to estimate

the arrival of the solute pulse from sparse data.

The program, CXFIT (Parker and van Genuchten, 1983) was used to

perform a nonlinear least-square fit of the analytical solution to the

solute concentration time series collected at each location. In this

program the retardation factor was set to 1.0 and the values of velocity

and dispersion coefficient were determined by minimizing the sum of the

squared residuals. Initial estimates: v 10 cm hr-~ and D = 50 (cm hr-1)2

were assigned to all locations. The coefficient of determination, r2,

which measures the degree to which the model fits the data, was used to

determine how well the selected model accounts for variation of the data.


The results of fitting the analytical solution to the individual

time series varied greatly (Figure 5). The coefficient of determination

ranged from .05 to .99. (The complete data set for Experiments I, II, and

III is presented in Tables 22, 23, and 24 respectively in APPENDIX C).

For each experiment there were locations where the analytical solution fit

well and other locations where the solution fit poorly. In the majority

of locations the estimate of pore water velocity coincided well with the

estimate of solute pulse velocity determined from the arrival time of the

peak solute concentration (Figures 5a and 5b). In approximately ten

percent of the sample sets double peaks occurred and the estimate of

solute pulse velocity became an average between the arrival times of each

peak (Figure 5c).

It is apparent that the model significantly underpredicted the

amount of solute transported in both Experiment I and III. The modeled

results for Experiments I and III could have been improved by estimation

of the indigenous nitrate and scaling the measured values. However, the

important aspect of this analysis was to improve the estimation solute

pulse velocity which was adequately achieved without scaling. For iodide

transport in Experiment II the modeled peaks were very similar to measured

peaks for a majority of sample sets (Figure 6). In the best case where r2

was greater than 0.95, such as at node 3-5 (Figure 6a), the model

accurately predicted the peak location and pulse shape. The description

was good within the uncertainty due to the sparse sample size. In the

average case, such as at node 2-7 (Figure 6b), the model predicted the

arrival of the peak concentration within one sampling interval.

Figure 5. Soil solution nitrate concentration time series and fitted model
from CXFIT with coefficient of determination (r2) for model fit for three
lysimeters (identified by E-W coordinate and N-S coordinate) for
Experiment I.

a) Location

b) Location

c) Location

14.25, 2: r2 0.96.

15.50, 2: r2 = 0.70.

13.50, 2: r2 0.35.

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

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