LEACHING MODELS FOR SUBSURFACE POLLUTION ASSESSMENT IN
AGROECOSYSTEMS
Mohamed M. Hantush1, Associate Member ASCE, Zhonglong Zhang2, Victor
Murray3, and Miguel A. Marino4, Member ASCE
Abstract
Unrestricted use of pesticides in agriculture threatens ground-water resources and
can have adverse ecological impact on the nation's receiving surface waters. In this paper,
we develop mass fraction models for exposure assessment and the regulation of
agricultural organic chemicals. The models are obtained by applying the method of
Laplace transformation to solute fate and transport equations. The models describe residual
mass emissions of pesticides below the root zone, to the water table, and in aquifers. They
emphasize the physical and biochemical processes responsible for the fate and transport of
organic solutes in the subsurface and their relationship to chemical properties and a set of
environmental factors. The potential of the use of the mass fraction models in the
management of agriculture pesticides is also investigated. A combined modeling and
ArcView GIS framework is used to assess soil and groundwater vulnerability to the
pesticide dicamba (commonly used for soybean and grains) in an agricultural watershed in
the Mid-Atlantic Coastal plain.
Introduction
Pesticides used in crop production are a major source of nonpoint source pollutants
to ground water, and their discharge to the nation's surface water may be a contributing
factor toward the decline of the living resources and the deterioration of ecosystems. Cost-
effective assessment tools are needed to regulate the use of agricultural chemicals, identify
areas which are potentially vulnerable to nonpoint-source pollution, and support ecosystem
restoration goals by improving the nation's water quality. Physically based environmental
simulation models can be cost-effective tools for resource managers when compared to
costly and prolonged field monitoring strategies. In this effort, we model solute leaching in
two distinct soil compartments, the soil-root layer and intermediate-vadose zone, and
implement them within ArcView GIS to agricultural watersheds in the Mid-Atlantic
coastal plain.
1 Hydrologist. Subsurface Protection and Remediation Division, National Risk Management Research
Laboratory, ORD, U.S. EPA, 919 Kerr Research Dr, Ada, OK 74820 (Corresponding Author)
2System Analyst, ManTech Environmental Technology, Inc., 919 Kerr Research Dr., Ada, OK 74820
3GIS Specialist, ManTech Environmental Technology, Inc., 919 Kerr Research Dr., Ada, OK 74820
4Professor, Department of Land, Air and Water Resources, and Department of Civil & Environmental
Engineering, University of California, Davis, CA 95616
1
Hantush

-------
Leaching models or indices are used for screening organic chemicals relative to
their mobility in the soil (Laskowski et al., 1982; Jury et al:, 1984Rao et al., 1985;
Mahmood and Sims, 1986; and Meeks and Dean, 1990). Beltman et al. (1995) developed
mass fraction models that describe leaching in the soil and subsequent mixing and
attenuation in ground water. Although index models are not suitable for predicting
concentrations, they nevertheless are simple and require data that are generally available
from scientific literature and soil survey reports. Further, they can be integrated with GIS
and produce an effective assessment tool for the management of nonpoint source in
agricultural watersheds (e.g., Khan and Liang, 1989, and Loague et al. 1995).
Leaching Below the Root Zone
Fig. 1 illustrates the conceptual soil ground-water system that is modeled here. We
assume soluble fraction of pesticide mass per unit area of soil, M0 [Kg], is instantaneously
mobilized at the surface by infiltrating water into the root zone. This may be represented
mathematically by a Dirac-delta pulse of input mass for which the average solute
concentration in the root zone is governed by the following equation (Hantush and
Marino, 1997), written in the differential form
dCr + prCr dt = -^-8(t) dt	(1)
QrhRr
where Pr=[1 + ln(2)(rr /A)(l + //)]/Tr; Cr(f) = average concentration in the root zone
[Kg/m3]; Tr= h R,J(v lQr) is the residence time in the root zone [day]; v* = percolation
below the root zone [m/day] - equals to precipitation plus irrigation rate minus
evapotranspiration,
Volat. vM„
ET
Rbot zone
Dscay
v Mr
Intcrviadosc soil
v Mu
Aquifer
Mg(x)
'///////////////////////////A
x=lx/2	X
*0
Fig.l. Conceptual soil-aquifer system
2
Hantush

-------
ET; |j, = (F S+c/h)X/[0.693 RrOr], Rr = liquid-phase partition coefficient, R,= 1 +(puA't/+Kr
Ku)l 0r\ a = KrKHDg/d; 0r = average volumetric water content in the root zone; h = depth of
the root zone [m]; pb = bulk soil density [Kg/m3]; X = pesticide half-life [day"1]; Kd =
distribution coefficient [m3/Kg]; Kr = volumetric air content in the root soil; KH=
dimensionless Henry constant; S = transpiration rate [day"1]; F = transpiration- stream
concentration factor; Dg = gaseous diffusion coefficient [m2/day]; and d= thickness of air
boundary layer on soil surface. Equation (2) is based on the following assumptions: i) well-
mixed root zone; ii) volatilization from soil surface occurs through an air boundary layer of
thickness d; iii) first-order rate reaction; and iv) passive plant uptake - rate of uptake is
proportional to dissolved solute concentrations.
One-dimensional transport and fate of solutes in the intermediate-vadose zone may
be described by the following boundary value problem
o,R,^r = ~-o,R,KC,	(2)
ot dz
Boundary conditions:
B.C. 1	Fu(0,t) = vCr	(3)
B.C. 2	C„(oo,0 = 0	(4)
Initial condition:	Cu(z,0) = 0	(5)
dC
where Fu is the solute flux, Fu = -QUDU —- + v Cu; Du= (Ku/Qu)KHDg+Dz is the effective
dz
liquid-phase diffusion coefficient [m2/day]; Dz = vertical dispersion in the soil [m2/day];
9U= volumetric water content in the intermediate-vadose zone; Rn= 1 +(puf(c&KuK[[)/ 9U; ku =
0,693/a is the decay rate coefficient [day"1]; and icu = volumetric air content in the
intermediate-vadose soil.
Leaching per unit area of soil, l(z,0, at depth z below the root zone at time t is
defined by:
Iu(z't) = iFu(z^)d^	(6)
The application of the Laplace transformation to Eqs. (1) and (2), and the use of Eqs. (3)-
(5) yield
~	~	i/l+4Tu P~\v+ku)-\ 1 Pu —
Fu (z; p) = Fu (0; p) e	H , and
v*/e 1	^
K(0;p)=mo -
hRr $r+p
~ r ~ t
where f{z\p)=\ f (z;t) e p dt, is the Laplace transform of the function fiz,t), Tu= H
JO
i?M/(v*/0u)is the residence time in the intermediate-vadose zone, and Pu = H (v* / 0U)I Du is
a Peclet number in the intermediate-vadose soil. The Laplace transformation of the
leached-fraction function Iu(z,t) follows immediately by applying the transform to (6):
Iu (z; p)= (1/ p) Fu (z; p) . Thus, residual solute mass that eventually leachs below depth z
3
Hantush

-------
is given by Iu (z) = lim 0 p Iu (z;p) - lim p^0 Fu (z; p), in which the substitution for
p^> 0
Fu (z; p) from (7) and taking the limit yields
'.(*) = ¦


1+4 kuTu p-
-up„-
l + ln(2)(rr/A,)(l + n)	(8)
The first of the products on the right-hand side, /r, accounts for the leached mass fraction
below the root zone.This equation describes the interactions among the physical and
biochemical processes and their net effect on vertical leaching in agricultural soils. Total
emissions of the pesticide to the water table, /„, can be obtained by the substitution ofH for
z in Eq. (8). Note that X accounts only for biodegradation only, the effective half-life ti/2,
which also accounts for the effect of roots uptake and volatilization from the soil surface,
is, thus, given by tm = A,/(l+[i).
The evaluation of Eq. (8) requires average values for each of dr and 0U. While 0r
may be assumed to attain the value at field capacity (e.g., Rao etal., 1985), average 9U can
be estimated assuming quasi-steady gravity flow: v = K(0U) and using appropriate
hydraulic conductivity-matric potential function (e.g., Clapp andHornberger, 1978). An
estimate for expected groundwater concentration may be given by Iu /(n B); in which n =
aquifer porosity, and 5 = thickness of the contaminant plume [m]. More details are given
later in the application section for estimates of groundwater concentrations on the basis of
multiple applications.
Emissions to ground water
Pesticides may undergo further attenuation along groundwater pathways because of
dispersion and biochemical decay in the aquifer porous matrix. The residual solute mass
convected across a plain normal to groundwater flow and below an agricultural plot of
dimensions 4 x ly (Fig. 1) can be obtained by the integration of Hantush and Marino
(1997) solution, which describes two-dimensional and depth-averaged concentrations in
ground water due to a rectangular nonpoint source, over the source area and in time:
M(x)I(MJJv) =
1
p*4y
V^+1

4r~ i
\-e
¦P (x/lx+1/2))
(9)
\-e
2
(x/lx-\l2))
x <1 / 2
where y = 1+ 4 Tg k/P", p = [(y1/2-l)/4] P"; P*= u IJDX is a Peclet number; Tg = lxRg/u is
ground-water residence time [T]; and u is the average pore-water ground-water velocity
[L/T], Eq. (9) describes the integrated effects of dispersion and biochemical reactions on
the residual mass flux of the contaminant past an aquifer section below and at a distance x
downgradient from the center of the agricultural plot. It can be used to estimate residual
mass of the pesticides entering a surface-water body (e.g., stream) downgradient from the
4
Hantush

-------
plot.
Implication on Management of Pesticides
A strategy for the management of organic compounds can be implemented by
imposing an upper limit F on the leached fraction in order to screen out those pesticides
that have the highest potential for ground-water contamination: Ir < F*Ma, and Ir is the
integrated mass flux (leaching) below the root zone. The substitution of the first of the
products on the right-hand-side of Eq. (8) - recall, represents leaching below the root layer
- for Ir in the above inequality and solving for T,/a yields
T
— > 1.44
X
f\-F^
F
1
1 + |J,
(10)
To illustrate the aggregated effect of volatilization and uptake by natural vegetative
cover, we consider the case \x = 1 (typical to the pesticide bromacil, which has
characteristics suitable for greater uptake by plants) and \x =10 (typical to the highly
volatile heptachlor). Then for a leaching factor F = 0.1, the substitution for |j, in (10) yields
Tr > 6.5 X and Tr > 1.2 X, respectively. Volatilization has much greater impact than crop
uptake on the management process by extending the pool of environmentally acceptable
pesticide; i.e., by expanding the list of relatively long-lived pesticides that can be used.
Thus, the use of process driven screening models may have implication on cost-effective
management of pesticides.
Nonpoint Source Assessment (GIS)
The study site (Locust Grove) is paired-watersheds of the Chester River drainage
basin. The Chester River constitutes the southern boundary of Kent County, Maryland, and
lies within the mid-Atlantic Coastal Plain along the northeastern shore of the Chesapeake
Bay (Fig. 2). The land in the Locust Grove site is predominantly agriculture, used to grow
mainly corn and soybeans in an annual rotation with winter wheat. Loamy soils, which
range from sandy loam to clay loam, dominate the landscape in that area, and the surficial
aquifer at the site consists of sand and gravel of fluvial origin. The relatively shallow water
table and high conductivity of the aquifer favor greater vulnerability of ground water and
stream baseflow to agricultural chemicals. The soil data is digitized from Kent County Soil
Survey maps and complemented by soil hydraulic properties, including organic carbon
fractions from available literature (e.g., Carseletal., 1988). Figure 3a classifies the soil in
the study area according to hydrological groups. Figures 3b and 3c are Arc View GIS
display of estimates of concentrations of the pesticide dicamba (commonly used for
Wheat-barley-alfalfa and Soybeans) in the intermediate vadose zone, Irl{duH), and ground
water, Iu/(nB), respectively. The estimates are based on one application of 2.5 lb/acre and
(February-April) climatic and hydrologic conditions.
5
Hantush

-------
Washington D.'

Kent
County
'
Fig. 2 Map of the study area (Locust Grove, Kent County, MD)
The above equations for the estimation of soil and groundwater concentrations are
programmed, using avenue scripts, and the results shown in Figs 3a and 3b are based on
running the scripts within the Arc View GIS project. Figure 4a displays the nonpoint source
module (NPSModeling) built into the Arc View project menu. Some of the required input
parameters are independent of soil features (or polygons in Arc View GIS), such as the
temperature, precipitation, Henry constant, initial mass application, ...etc., and they are
entered for processing by the avenue scripts through a Visual-Basic-type user-interface
features as shown in Fig. 4b.
The soils in Chesterville Branch appear to be primarily class B, whereas they are
primarily class C in Morgan Creek (Fig. 3a). Thus, the former basin is relatively more
drained than the latter, thereby resulting in greater expected maximum concentrations of
dicamba in the intermediate-vadose zone and in ground water in the Chesterville Branch.
The effect of the spatial variability of soil texture and drainage properties, and the depths
to water table, is evident on the distribution of estimated concentrations in the
intermediate-vadose zone and in ground water. The excessively estimated concentrations
in the intermediate-vadose zone, especially toward the Sassafras River and in other
isolated areas, are attributed to reference surface elevations where the depths to the water
table approach zero, H= 0 m (Fig. 3b). Errors in kriged depths to the water table may also
be a contributing factor. We emphasize that the estimated concentrations in the soil and
ground water only present expected magnitudes and are short lived, especially those
estimated in the soil (Fig. 3b), since they are associated with a single application of the
pesticide. However, the estimated concentrations in ground water may persist for a longer
period, especially for multiple applications of dicamba. If a pesticide is applied at a certain
frequency, to (in this application to = 1), a quasi-steady state may develop and the
concentration in ground water can be estimated from
6
Hantush

-------
(a)
HYDROLOGICAl SOIL GROUPS
|yliBWE&
i
j^Mlilni.t;,
* #Wf
¦iilasi^^^«a»
¦i
mm/mm
^¦¦Mp
¦L,
lllfiJB|l
WSm

JBHl
ishksk
mm
KATES Of
JNBlTfWfON
«lliiiiis
¥ ;,r
A
(b)
i
Morgan Creek
£».
%
¦A
-\
;}
¦ • 5? 55 57 55-n?0* 137 04. 513«*)
0.W • M» 1 TS-JZfi
5«- u 17 U.17HJ5.M
Fig. 3 (a) Soil hydrological classes; (b) Estimated concentrations (ppb) of dicamba in the
intermediate-vadose zone in Spring; (c) Estimated concentrations (ppb) of dicamba
in ground water in Spring (the intermediate-vadose zone in Spring
C\j) =
®IU (j)
n B
(11)
in which C (j) is an estimate of quasi-steady concentration in ground water (Kg/m ) and j
is an index, which accounts for seasonally-averaged hydrologic and climatic conditions.
An alternative approach for calculating the quasi-steady state groundwater concentrations
is to assume that applications of a pesticide are flushed into the water table during a given
season by the cumulative recharge (in volume) (Beltman et al., 1995):
=	(12)
P(j)
7
Hantush

-------
in which P(j) = is the cumulative recharge during season j [m3].
As Fig. 3c shows, biodegradation of dicamba in the intermediate vadose zone
(assuming half-life = 50 days) and mixing in the aquifer (B = 5m) have the potential of
significantly reducing (by orders of magnitude) dicamba concentrations in ground water.
It is often the case that half-live values (X) of pesticides that are reported in the
literature lack sufficient details regarding the conditions under which they are measured,
therefore, are subject to uncertainty. Information as to whether measured values of the
half-lives were made in the laboratory or under filed conditions, is rarely sited in the
literature. Field-measured half-life usually incorporate the combined effect of all possible
loss pathways (biodegradation, volatilization, and crop uptake), and their use in physically
based resolute models is bound to overestimate actual losses, and consequently,
underestimate the magnitude of potential pollution problems. In case of uncertainty in
reported values of X, the effective half-life ti/2 = X/( l+|_i) can be treated as a random
variable. And losses due biodegradation, volatilization, and crop uptake, are lumped into
one parameter, ti/2, rather than being evaluated individually (recall, data requirements for
the estimation of S, a, and |_i).
Summary and Conclusions
A methodology is presented for the assessment of groundwater vulnerability to
pesticides in agricultural watersheds. The methodology integrates simple but physically
based analytical models with ArcView GIS to assess the fate of pesticides in the
subsurface. The models describe the combined effect of the physical, (bio)chemical
processes, and related environmental factors on the fate and transport of the organic solutes
in the soil and ground water. The models are capable of describing leached mass fractions
of solutes below the root zone, and the residual mass emissions to ground water. The
leached mass fraction and residual mass emissions are estimated on the basis of cumulative
mass flux past a particular section normal to the flow in the subsurface porous media. The
resolution of the models is such that processes of crop uptake and volatilization of the
pesticides can also be accounted for as well as advection, dispersion, equilibrium
adsorption, and biodegradation. It is demonstrated that the ability to quantify potential loss
pathways by crop uptake and volatilization may leads to improved (more cost effective)
management of the pesticides, with the objective of protecting groundwater quality.
The modeling-GIS framework is applied to mid-Atlantic coastal plain agricultural
watersheds for the assessment of potential impact of the use of the pesticide dicamba on
groundwater quality. The primary objective of using ArcView GIS is to investigate the
impact of the spatial variability of soil hydraulic, biochemical properties and other
environmental factors on the variability of dicamba emissions to ground water. Preliminary
estimates indicated that ground waters in the relatively well-drained areas around
Chesterville Branch received higher doses of dicamba, probably due to relatively shorter
residence time in the soil where (bio)chemical degradation has more impact than in the
relatively poorly-drained areas in the Morgan Creek watershed. Further assessment is
needed because of uncertainty in soil and hydrologic parameters and estimation errors in
kriged depths to the water table. The input parameters were based on measured values and
8
Hantush

-------
average values reported in the literature.
Estimation of the impact of uncertainty of input parameters on the variance of the
estimates of concentrations is currently under investigation, using first-order and second-
order analysis.
Notice: This paper has been reviewed in accordance with the U.S. Environmental
Protection Agency's peer and administrative review policies and approved for presentation
and publication.
References
Beltman, W.H.J., J.J.T.I. Boesten, S. E. A. T. M. van der Zee. 1995. "Analytical modeling
of pesticide transport from the soil surface to a drinking water well." J. Hydrol. 169,
209-228.
Carsel, R. F., R. S. Parrish, R. L. Jones, J. L. Hansen, and R. L. Lamb. 1988.
"Characterizing the uncertainty of pesticide leaching in agricultural soils." J. Contam.
Hydrol. 2, 111-124
Clapp, R. B., and G. M. Hornberger. 1978. "Empirical Equations for Some Soil Hydraulic
Properties." Water Resour. Res. 14(4), 601-604.
Hantush, M.M., and M.A. Marino. 1997. "An analytical model for the assessment of
pesticides exposure levels in soils and groundwater." Journal Environmental Modeling
&_Assessment, 1(4), 263-276.
Khan, M.A., and T. Liang. 1989. "Mapping pesticide contamination potential.
Environmental Management 13:233 -242.
Jury, W.A., W.J. Farmer, and W.F. Spencer. 1984. "Behavior assessment model for trace
organics in soil: II. Chemical classification and parameter sensitivity." SoilSci. Soc. Am.
Proc., 13, 567-572.
Laskowski, D. A., C. A. I. Goring, P. J. McCall, and R. L. Swann. 1982. Terrestrial
Environment, p. 198-240. In: R. A. Conway (ed) Environmental Risk Analysis for
Chemicals. Van Nostrand Reinhold Co., NY.
Loage, K., R. L. Bernknopf, T.W. Giambelluca, and R.E. Green. 1995. "The impact of data
uncertainty upon regional scale leaching assessments of non-point source pollutants."
Applications of GIS to the modeling of non-point source pollutants in the vadose zone,
ASA-CSSA-SSSA Bouyoucos Conference, Mission In, Riverside, CA.
Mahmood, R. J. and Sims, R. C. 1986. " Mobility of organics in land treatment systems."
J. Environ. Eng., 112(2), 236-245.
Meeks, Y. J. and J. D. Dean. 1990. "Evaluating Ground-water vulnerability to pesticides."
Journal Water Resources Planning andMangement. 116(5), 693-707.
Rao, P.S.C., A.G. Hornsby, and R.E. Jessup. 1985. "Indices for ranking the potential for
pesticide contamination of groundwater." Proc. Soil Crop Sci. Soc. Fla., 44, 1-8.
9
Hantush

-------
(a)
AirVii;w lir; V).isiiin "I Ita
Eiie £jit	Them* Jjraphici ^ndow
a b wem\ a:
«>-
IIJ'I'.M
Tables:'' ;
'L
€>
Layouts
fgii
HlH
.f
Hal
Help
Complete Mixing Calculation
Mass. FracbOfiRooi Zo^el
Mass Emissio.nlVados-e!
. •••-" ^
Mass F:action(RoofZone|
JL
ffli *1
-tie 1:|
JSIllf:
vjuii i ypca
(b)
1 r«r.u\t 1j|i»v»- J'|m|i'. {
7 e-'ipreture TIGC} 15 3-1
Henij,- constant: khstar . f jJALOJ'P
Gaseous diffusion coe!: Dg!c:2'di |C.«32
0.002
- ] 1	
¦R^feidsifal-te: i'amda[daff f50
r.."" *		
j
FiaciOs'-aOofs: FWseason] i 1J 36
^p3tran^ipt!Qri^f^(iri^a^:|2G2
I nniil i jlil>>; Phi|.:(J
E.'tfer ail xer defined data. click CK to continue-
Groundwater piume fhiokriess' E[aj 15
Leachhg facio? Fl i i 25
T inspiration reduction factor: gar:-"s jCf
Lea!areaircie.i ; jCi
initial r,i3?s:M0(kg/rs:i 10 CC.J2SC155
Aqjife' pGrosit.y 1 0 i
Potential Evapcti: ETu%'dayJ C HT1
MCL HalipFbi m
'•¦"e-'jcai dispeitlv'lo: AbhaLlcm] j:
Fig. 4. Graphical user-interface features for the Nonpoint Source Assessment Module in
Arc View GIS
10
Hantush

-------
Stratified Soil Profile
The semi-infinite solution of (8) is based on the assumption that dispersive flux is
negligible when compared to advection. Shallow impeading layers (e.g., silty-clay layer or
loess) can slow down convected solute fronts and produce a significant upward diffusive
flux. A semiinfinite solution of Eq. (2) is not valid in this case.
The first equation in (7) is a general Laplace-transform representation of leaching
in a distinct soil layer in a stratified soil. Starting with the upper most layer and assuming
that a semiinfinte solution is valid, one can deduce that leaching below a stratified soil
profile made of N layers is given by
In
f I ^
yMoJ
1
N
^ 7=1
l + 41„(2)i-(l + n()p;
-1
p.:
(9)
in which, Tu' and PL1' are the apparent residence time and the Peclet number of layer i,
respectively, and \Xi = 0 in layers where losses by root uptake and volatilization are
negligible. This equation may be used to described leached residual mass fractions of
pesticides in stratified soils. It accounts for the heterogeneity of the soil profile and spatial
variability of its hydraulic and (bio)chemical properties. An estimate for expected
groundwater concentration may be given by I\/(n B); in which n = aquifer porosity, and B
= thickness of the contaminant plume [m].
11
Hantush

-------