SPATIAL VARIABILITY OF SOLUTE LEACHING IN A SANDY SOIL
By
ERIC GEORGE FLAIG
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1990
ACKNOWLEDGEMENTS
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.
TABLE OF CONTENTS
page
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 ShortRange and
LongRange Variability ......................................... 177
APPENDIX A NITRATE AND IODIDE ANALYSIS USING ION CHROMATOGRAPH.... 179
APPENDIX B EVALUATION OF TENSION LYSIMETERS ....................... 182
APPENDIX C COEFFICIENTS FOR REDISTRIBUTION MODEL.................. 198
LITERATURE CITED .................................................... 205
BIOGRAPHICAL SKETCH ................................................. 225
LIST OF TABLES
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 twodimension 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 hr1) for three field Experiments: I, short
transect; II, grid; III, transect............................. 118
9 Statistical properties for estimated solute mass loss (SML,
kg ha1) 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 (SMLres) 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 20mg 11 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 20mg/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
vii
LIST OF FIGURES
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 EW coordinate
and NS 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 EW coordinate
and NS 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
viii
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 B85, b) site B82.. 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
(SMLresid) 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
(SMLresid) 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
2dimensional 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 2dimensional random field:
a) mean averageerror, b) variance of averageerror........... 164
26 Mean squarederror for 50 sample sets of increasing sample
size for different patterns taken from a 2dimensional random
field: a) mean mean squarederror, b) variance of mean
squarederror................................................. 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
SPATIAL VARIABILITY OF SOLUTE LEACHING IN A SANDY SOIL
By
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 surfaceapplied solute
pulses (SPV) were determined by fitting an analytical solution to the
convectivedispersive 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 09%
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 twodimensional random field with known covariance was simulated
using turningbands 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 crossvalidation 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.
xii
INTRODUCTION
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,
1983).
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, highvolume
irrigation representing design storm events likely to result in potential
2
offsite 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 steadyunsteady and saturatedunsaturated 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 microrelief 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 highvolume 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 onedimensional transport model was to estimate
4
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
determined.
REVIEW OF LITERATURE
Introduction
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 surfaceapplied solutes during infiltration
can be estimated by assuming pistonlike displacement of initial water by
the invading water (Warrick et al., 1971; Kirda et al., 1973; Ghuman et
6
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
depthdependent 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 bypassed by the applied water. In either case,
water applied following the solute pulse determined the leaching depth.
7
For longterm 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 timematched infiltrationredistribution. 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
8
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.
ConvectiveDispersive 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 convectivedispersive 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 11), with time
(t) and depth (z) where water content 9 (cm3 cm3), pore water velocity v
(cm hr1), and the dispersion coefficient D (cm2 hr1) 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
9
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).
10
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 80fold 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
11
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 convectivedispersive 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
12
term transient flow. Under bare soil with essentially no immobile water
the pistonflow 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
below.
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 semiinfinite 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
13
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,
1982).
Analytical solutions have been fit to both residentaveraged and
fluxaveraged solute concentration time series. Fluxaveraged solute
concentration samples are those collected from actively flowing water.
These samples are typically collected from the discharge of soil columns.
Residentaveraged samples provide a measure of solute concentration in all
water resident in the soil. Residentaveraged samples may be collected by
pore water extraction techniques.
An analytical solution for CDF for fluxaveraged 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)
az
S o 0 < t < to
cf(O,t) = (4)
0 t > to
The boundary condition (4) represents a flux boundary, a thirdtype 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,tto) t > to
where
A(z,t) erf Rzvt + t exp(+i') erfc Rzvt (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 infiltrationredistribution (Wagenet et al.,
1976).
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
15
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
16
of the 9hK 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
17
used successfully to describe nitrate leaching from relatively uniform,
unstructured soil without preferredflow 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 / (e8e1), (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
18
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, (Ofefc) 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; AmoozegarFard 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
19
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 lognormally
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
20
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 lognormally 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,
AmoozegarFard 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 firstorder approximation to
predict the effect of spatial variability on water flux. This method
21
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)
0
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
22
CIN with constant water flux, the solute concentration distribution across
the field at L is as follows:
CL(I) = CIN(II') fL(I') dl' (10)
0
The inlet concentration arrives at Z=L unchanged because dispersion is not
considered. Dagan and Bresler (1979) explicitly included dispersion by
replacing CIN(II') with the appropriate solution to the onedimension 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(II') 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.
23
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.
24
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 onedimensional systems assuming all soil properties
except K(O) remain constant and Ks was randomly selected from a lognormal
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 infiltrationredistribution (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
25
simplified method for determining soil water flux during infiltration and
redistribution could avoid this difficulty.
Infiltration
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.
Redistribution
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
redistribution.
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
27
defined a "shortterm" 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 halfhour 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
28
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.
Soilwater 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 coarsetextured soils than
for finetextured soils (Hillel, 1980). Staple (1966) reported that water
content during redistribution in a sandy soil was approximately constant
a) soil water flux
20.0
15.0
10.0
5.0
I I t t I I I I I I
1 2 3 4
Duration, t
5 6 7 8
(hr)
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)
L.
..C
E
U
X"
_D
L
0
0)
Eustis sand
.716
 q = 2.889 t
 
_1
__ _
30
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 massbalance or estimated by using Richards
equation. Solutions of Richards equation require estimation of the soil
moisture characteristic and soil water conductivitypressure head, K(h).
Estimation of these data requires considerable field monitoring.
Soil water flux through the profile can be determined by calculating
a water massbalance 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
31
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 soilmoisture 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
32
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
33
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 longterm 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 zerotension
lysimeters are the most common methods for soil solution sampling. The
zerotension 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. Zerotension 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 poresequences 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 nonadsorbed 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
35
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 AmoozegarFard,
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 fieldwide frequency distributions best described by
the lognormal 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
36
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 5060% 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 lognormal
37
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,
38
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, soilwater
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
1520 m and saturated hydraulic conductivity (Ks) a range of 1520 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 longterm 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 microrelief.
a) components
2.0 
1.5
1.0 
nugget,
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
Spherical
  Gaussian
 Exponential
 Linear
Lag Distance, h
Figure 2.
models.
Theoretical semivariogram: a) components, b) semivariogram
variance sill,c
model, \A
range,a)
C
40
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 soilwater 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.
41
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 s1. 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).
42
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
43
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 soilwater
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
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
44
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 kth 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(x2m(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
45
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 secondorder 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
realizations.
Semivariogram
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.
46
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 Mestimator (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
47
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 crossvalidation procedure
(Gambolati and Volpi, 1979). Crossvalidation 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 nl values using the kriging procedure. This is repeated for
each sample point providing nl 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.
48
The accuracy of the model in a meansquared 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:
(23)
RMSE 1/n e2 1
i1 
Oji
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). Crossvalidation 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 crossvalidation.
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
the
low
model
kriging
smoothly follows the experimental
error (Yost et al., 1982b). Kriging
49
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 nonzero value, kriging error will be significantly
understated.
High crossvalidation 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).
Kriging
The third step to the geostatistical process is kriging. Kriging is
a linear, unbiased estimation process for estimating point or block values
50
(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:
n
Z*(xo) = Z Ai Z(xi) (25)
i1
that is unbiased if
E[Z*(xo)] E[Z(xo)] (26)
which requires the following condition on the weights,A,
n
E Ai = 1 (27)
i=1
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(xixj) (28)
i i J
where: y(xi,V) average semivariogram between points and domain V
(V,V) average semivariogram of domain V
y(xixj) 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":
n
Zi y(xj,xi) + = (xv,xi) i= 1,...,n (30)
i=1
n
and subject to: Z Ai 0 (31)
i=1
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:
n
aoK = Z Ai 7(Xj,Xi) + p (32)
i=1
52
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 (xixj) (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(xixj) semivariogram among all points.
The kriging variance for block kriging is given in equation (34).
a= Aik iu + f uu(34)
i
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 withinblock variance is removed from the
error term.
53
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 nonnormal (Rendu, 1980). These alternative
methods provide useful techniques in limited circumstances.
Drift
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 largescale, 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,
54
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 longrange 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 29point moving average.
Hamlett et al. (1986) removed drift from soilwater 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
55
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 firstorder autoregression model to describe soil
temperature. They used a fourthorder 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 meansquared prediction error using cross
validation.
Where the drift function or the semivariogram are complicated, these
functions may be simultaneously identified by utilizing intrinsic random
56
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 kthorder IRF is defined
as a random process that requires a korder 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 IRFk (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
57
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 RodriguezIturbe (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 wellbehaved phenomena with CV less than one, kriging
58
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,
1981).
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
59
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
properties.
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.
Uncertainty
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
60
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 firstorder 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.
61
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 secondorder 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
analysis:
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 firstorder 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 datacollectionnetwork design has the objective of
altering the network configuration until minimum error is achieved (Bras
and RodriguezIturbe, 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.
63
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
datacollection 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
64
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
spatiallydependent 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 9fold 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
65
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 userdetermined
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
representation.
Hughes and Lettenmaier (1981) used an iterative algorithm to adjust
the location of sampling sites in order to minimize regional kriging
66
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 twodimensional 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.
Simulations
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 twodimensional field on known covariance can
be generated. This can be accomplished using the turningbands 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 twodimensional 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 twodimensional 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).
Synthesis
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,
68
extrinsic sources, particularly variability in application rate may be
important.
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 convectivedispersive flow model for transient
flow.
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.
69
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 nonconstant 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 crossvalidation is the
best available technique. Rather than optimizing the model fit to all
data, crossvalidation 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, intersample
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.
METHODS AND MATERIALS
Introduction
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
71
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 convectivedispersive flow (CDF) equation was
fit to the time series of solute concentration from each location using a
nonlinear leastsquares 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.
72
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.
Laboratory
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 ceramiccup
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
74
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
design.
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 thinwalled 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,
10cm depth increments. From these samples the following analyses were
conducted: bulk density by core method (Blake and Hartege, 1986) (Table
1); particlesize 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 cm3
Ave. SE
5.20
4.67
5.28
5.05
5.46
0.22
0.11
0.06
0.10
0.12
15.1 0.32
14.2 0.25
1.65
1.60
1.54
1.52
1.49
0.036
0.030
0.021
0.016
0.010
Table 2. Soil physical properties of typical Kendrick sand (Loamy,
siliceous, hyperthermic, Arenic, Paleudult).
Particle Size Distribution (% < 2mm)
Depth
(cm)
VC
(dia. 2.0
mm) 1.0
C
1.0
0.5
M
0.5
0.25
F
0.25
0.10
VF
0.10
0.05
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 cm3)
(SE)
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
Depth
 cm 
12.4
0.11
Field Leaching Experiments
Three experiments were conducted during 19851986 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 intersample 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 highvolume 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
77
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
redistribution.
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 2hour, 10year
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 25ASFPTNT
with 0.32 cm nozzle) located along the sides of the plots outside a 2m
buffer strip. The sprinklers were clustered in groups of four on top of
fivefoot risers. The water supply was received through a 5cm 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 16m 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
S8o
1 8
meters
o o o
0 0 O
***o**:***
o
0
** ** *
For Detail See Inset
* o Lysimeters, 1.0 m Spacing
b) Areal view of lysimeter placement for Experiment I.
Inset
0 meters
0 0 0 0 0
0ooooooooooo o olololololololololo o 0 o
0 0 0 0 0
8 20
Legend
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.
ooD
o 0 0
o 0 o
0 0 0
o 0 0
0 o o
0 0 0
0 0 0
I
c) Areal of lysimeter placement for Experiment II
Experiment II, Square Grid
0r4
meters
o Vertical Lysimeters at 1.0 m Spacing
d) Areal of lysimeter placement for Experiment III.
Experiment III, Long Transect
0
0
0000000 000
0000000000
meters
Figure 3. Continued
5.
09 0 0 0
0
.0 0 0 0 0
10
0
000
0
0
0 0 0
0 0
0
0
0 0 0
0
0o
000
80
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 250ml opentop 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 ha1.
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 5cm long and 0.2cm 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
Duration
Chemigation hr 0.62 0.33 0.33
Duration
PostChemigation hr 2.34 3.91 2.83
Irrigation
Average cm hr1 4.60 4.61 3.58
Irrigation Rate
Solute Applied Nitrate Iodide Nitrate
Average Areal kg ha2 5.61 12.2 5.15
Chemigation Rate
Solute Concentration mg 11 17.7 80.3 43.6
in Chemigation
Number of samples 17 12 20
Collected per site
with nominal 1bar bubbling pressure (Soil Moisture Equipment, Inc.) to an
open end of a 70cm length of 0.20cm 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.05cm 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).
82
Each sampling manifold consisted of a 50cm 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 50ml 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 nonreacting
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 2m 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
Vacuum
To Pump
A
1/4 Nylon
Tube
3 5m x .05cm Nylon Tube
. (9) Additional
Collection Tubes
50ml
Polycarbonate
Tube
Sampling manifold and collection tubes for tension lysimeters.
Pinch
Clamp
Tension
Lysimeter
Porous
Ceramic
Cup
o:n, 11 4
'
.:.
:' : .: '

r

Figure 4.
84
using a thinwalled 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
samples.
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 reapplied. This resulted in an
initial purge of water from the cups which was discarded. Vacuum was re
85
applied for approximately 1015 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
vacuum.
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
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.05cm 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 lagtime 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 fluxaveraged 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 leastsquare 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 hr1)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.
87
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 35 (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 27 (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 EW coordinate and NS 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.
