UNCERTAINTY IN LEACHING POTENTIAL OF NONPOINT SOURCE POLLUTANTS WITH
APPLICATION TO A CIS
Mohamed Hantush
Hydrologist
U.S. EPA, National Risk Management Research Laboratory
Zhonglong Zhang
System Analyst
ManTech Environmental Technology, Inc.
Miguel A. Marino
Professor
University of California, Davis
Enamul Hoque
Scientist
ManTech Environmental Technology, Inc.
ABSTRACT
This paper presents a stochastic framework for the assessment of groundwater pollution potential of nonpoint source
pesticides. A conceptual relationship is presented that relates seasonally averaged groundwater recharge to soil
properties and depths to the water table. The analytical relationship shows a linear association with the soil saturated
hydraulic conductivity, and predicts less recharge for shallower water tables. A modified index model for pesticides
mass emissions in the soil is developed, which accounts explicitly for biochemical degradation in the mobile and
immobile soil-water zones and losses by volatilization and crop uptake. The index model shows that the effect of
nonequilibrium transport on the residual mass emissions is dependent on ratio of the mass transfer coefficient and
the first-order degradation rate in the immobile phase. The stochastic framework utilizes first-order approximations
of the mean and variance of each of the recharge and the residual mass emissions in the soil. Soil and chemical
properties and related environmental factors, which affect the fate and transport of pesticides, and the depth to the
water table are modeled as random variables. The environmental-fate models are integrated with a GIS, and the
stochastic framework is applied to assess potential nonpoint-source vulnerability of shallow groundwater to the
pesticide dicamba in Mid-Atlantic coastal plain agricultural watersheds. It is shown that recharge estimates on the
basis of the SCS abstraction method resulted in lower expected dicamba concentrations in groundwater, than when
leaching is based on the conceptual model. This may be attributed to the fact that in the former leaching is averaged
over the soil hydrologic groups, whereas in the latter leaching reflects the spatial variability of the soil types and
environmental factors, including the depth to the water table. In the analysis, biochemical degradation of dicamba is
ignored below the root zone and the vulnerability results were based on a single application of 2.51b/acre. The first-
order estimates of the variance of dicamba mass emissions and estimated groundwater concentrations, were of order
of magnitude of their respective means and greater. While some of the mean groundwater concentrations showed
significant residual levels of dicamba, these values however should be viewed as subject to high uncertainty. Given
the expected uncertainties in the input data and model errors, regulatory decisions and environmental land-use
planning should take into account estimates of the uncertainties (e.g., variances) associated with predictions of
groundwater vulnerabilities.
INTRODUCTION
Pesticides are sited among the most common nonpoint source pollutants (NPS) pollutants, which threaten
sustainable agriculture and constitute a great threat to surface and subsurface drinking water resources (Loague, et
al., 1998). There is a heightened concern among the public and policy makers that soil and groundwater resources
are continuously threatened and degraded by improper use and mismanagement of agrichemicals. Proper land-use
-------
planning can circumvent potential harmful environmental impacts, such as groundwater contamination. Conceptual
environmental-fate models offer a promising avenue for the management of agrochemicals by forecasting the
movement of NFS pollutants and identifying potential loss pathways in soil, with the objective of preventing
pollution potential beforehand.
Modeling NFS pollutants, however, is a complex and a demanding task. The volume of information required for
temporal and spatial characterization of the parameters and variables needed, in even the simplest conceptual models
of solute transport in the vadose zone, is tremendous (Crowin et al. 1997). A Geographic Information System (GIS)
can be used as an effective tool for the management of volumes of data of spatial and temporal nature, and when
linked with environmental-fate models would make spatial data available for environmental analysis (Hoogeweg
and Hornsby, 1995). It is a tool designed to acquire, store, analyze, and display spatial information. The display of
spatially related information by a GIS can be very useful, especially in identifying spatial patterns of significant
values to regulatory agencies and policy makers for planning purposes. Khan and Liang (1989) were among the first
who applied GIS to create groundwater contamination likelihood maps for Oahu, Hawaii. They used the attenuation
factor AF of Rao et al. (1985) and Mahmood and Sims (1986) to assess the relative likelihood of groundwater
contamination. Rather than using complex mathematical models that describe the transport of pesticides in
subsurface systems (e.g., Wagenet and Huston, 1986; and Carsel et al., 1984), the simpler AF models (commonly
referred to as indices) are more favorable due to lack of independent measurements for parameters required in the
complex models. Kleveno et al. (1992) evaluated the attenuation factor (AF) by comparing it with the more rigorous
and conceptual pesticide root zone model PRZM (Carsel et al., 1984). They indicated that AF compared well to
PRZM in ranking the relative mobility of a selection of pesticides, in heterogeneous soil and under conditions of
variable recharge. Hoogeweg and Hornsby, (1995) created a graphical user interface for an integrated
GIS/environmental-fate models to facilitate access to multiple fate and transport models and the evaluation of
pesticide effects on groundwater quality. Loague et al. (1995) suggested the use of mobility indices as physical
factors in a regional integrated assessment approach (RIA) using GIS.
This paper emphasizes four particular aspects of groundwater vulnerability assessment including application to a
case study using a GIS. First, since recharge is the only mechanism responsible for the transport of pollutants from
the source to the receptor, its importance is recognized and a new conceptual relationship is presented for the
estimation of local temporally averaged (quasi-steady) groundwater recharge. The relationship is based on soil-water
retention model of Campbell (1974) and requires the knowledge of depth to the water table. Secondly, mass fraction
emission models are presented, which describe fate and transport of pesticides in soil, including the effect of
nonequilibrium transport due to mobile-immobile solute exchange by diffusion. It is worth noting that the utility of
existing indices (e.g., Rao et al., 1985; and Mahmood and Sims, 1986) is based on field based pesticides half-life
values, and their use under different field and climatic conditions may therefore be subject to uncertainty. Since
organic compounds, such as pesticides, differ in their characteristics and their behavior in the environment, then the
ability to quantify the relationship between their half-lives and the different loss pathways may reduce uncertainties
and improve model predictions. The relationship between effective (field scale) half-life and biochemical
degradation, volatilization, crop uptake, and losses in immobile water regions will be addressed. Thirdly, a
stochastic framework, using first-order analysis, will be developed for the estimation of means and variances of the
recharge and mass fraction emissions below the root zone, in terms of uncertainties in model input parameters.
Assessment models, such as the AF, are based upon soil, chemical, hydrologic, and environmental data that are
usually sparse and highly uncertain. Therefore, the impact of uncertainty of data, including errors in the conceptual
model, should be quantified as a measure of the reliability of predicted assessment maps that are based on fixed
input parameters (e.g., Loague et al., 1995). The methodology will be applied to assess shallow groundwater
vulnerability to the pesticide dicamba in paired agricultural watersheds in the Mid-Atlantic coastal plain. The
objective of the stochastic framework is to establish the extent to which the mean dicamba mass emissions are
reliable predictions of actual assessment maps of groundwater vulnerability. GIS displays of an environmental-fate
model outputs should be viewed with uncertainty and therefore has no real utility without a corresponding map of
associated uncertainties (Loague et al., 1998).
-------
DETERMINISTIC ANALYSIS
Conceptual Model For Quasi-Steady Recharge
In this section we develop an equation, which describes quasi-steady percolation below the root zone, as a function
of average depth to the water table and uniform (depth-averaged) soil properties. The Buckingham-Darcy flux law
can describe the flow in the unsaturated zone (Bear, 1972):
oz
(1)
in which, q(z,t) is the flux per unit area [m/d]; K(y) is unsaturated hydraulic conductivity [m/d]; \\i is the suction
head (negative of matric potential)[m]; z is the depth below the surface [m](Fig. 1); and t denotes time [d]. Note that
in (1) q follows the sign convention of assuming a positive value when it is defined downward.
We first define the time-average operator <> as
1 T
i
(2)
in which f(t) is some function of time, and T is a time interval over which (1) will be averaged and quasi-steady q is
defined. If we fix t at an arbitrary value t = T, 0 < T < T, and replace the partial derivative in (1) with a total
derivative (i.e., denoting derivative only with respect to z at fixed t), then Eq. (1) can be written as
dz
1
(3)
in which \\IT = \j/(z; t=x) and qx = q(z; t=r) are functions of z for a fixed t.
Precipitation.
ET
Runoff
Soil Surface
Fig. 1 Illustration of water balance in unsaturated zone.
Because abstractions due to roots uptake and evaporation are assumed to occur only in the root zone, q can be
assumed to be independent of z. The integration of Eq. (3) from the soil surface to the water table, z = h to z = H(T),
yields
-------
7!+ff(T) 0
j dz = j -
h vy,(T) I
where H is the depth to the water table below the root zone at t = T [m], and yr is the suction head at the interface
between the root zone and the intermediate vadose zone at t = T [m].
The relationship between the hydraulic conductivity and suction head for different soils may be describe by the
following equation (Campbell, 1974)
K(\v)=K\ , a= (5)
IvJ 26 + 3
where Ks is the saturated-soil hydraulic conductivity [cm/day]; \\is = saturation suction [cm]; and b is an empirical
parameter. Equation (5) is based on the following moisture characteristic
(6)
in which Q is the soil volumetric water content and 6S is the saturated water content. A data set for the parameters Ks,
0S, and b for different soils can be obtained from Clapp and Hornberger (1978). It should be noted that Eq. (6) is not
accurate near saturation (Campbell, 1974). The latter authors proposed a short parabolic section near saturation to
represent gradual air entry. The substitution of right-hand side of (5) for K in (4) and the application of the time-
average operator <> to the resulting integral equation yields
(7)
Where we have droped the subscript T notation from \\i in the integrand because \y is a dummy variable, and replaced
qx with q(x) to denote a functional relationship with respect to T. To the first order, the approximate integral equation
which governs quasi-steady recharge to the water table is given by:
o -
(8)
in which H = , yr = , and q = . Equation (8) is a first-order approximation to (7), since the left-
hand side of the latter is expressed in terms of time-averaged variables. The truncated higher order terms are
expected to be significant for large temporal variations in H(T), \J/(T), and q(x), and can be quantified using second-
order analysis in terms of the variances of fluctuations of these variables (second moments). However, this is beyond
the scope of this paper. Thus, Eq. (8) will be used to develop an approximate closed-form solution for quasi-steady q
in a uniform soil profile. After a simple transformation and expanding the integrand in (8) using Taylor series, it can
be shown that
(9a)
^VMM ^V ZCl+i Vr )
If the following conditions are satisfied
and H>\iir>\iis (9b)
-------
Its worth noting that this equation is consistent with the solution of the Buckingham-Darcy equation at equilibrium
conditions, since both predict that y(0) = \\ir = H when q=0. The condition \\is< v|/r is necessary because Eq. (5)
requires that \\is< \\i for K(\|/)< Ks to be a valid definition of unsaturated conductivity. Equation (9) is valid when the
water table has risen in response to a recharge process; therefore, its not valid for the specific case where a
decreasing recharge rate causes a declining water table. This equation predicts a linear association with the saturated
hydraulic conductivity and greater recharge for deeper water table, because greater recharge is required to elevate
deeper water tables. An interesting point may be that (9) provides an indirect mean for the estimation of percolation
below the root zone. That is, it rests upon the response of the water table to a net recharge process rather than flow
balance at the soil surface and the root zone. Thus, it may be used in case of uncertainty in irrigation practices and
actual evapotranspirative losses.
Leaching of Pesticides in the Soil
In this section we present a physically based relationship for leached mass fractions below the root zone and residual
emissions to the water table. The relationship considers the effects of biochemical degradation, dispersion, and
mobile-immobile nonequilibrium solute transport exclusively. Figure 2 illustrates the system under consideration.
Leached mass fraction of a pesticide mass M0 incorporated in a thin surficial soil layer may be described by the
following equation (Hantush et al., 1998)
(10)
where
R (rt I Is *m \ ^ m
(H)
ln(2)RrmQrm
where a = Krm KHD^d is the volatilization rate coefficient [m/d]; Krm= volumetric air content in the mobile phase;
KH= dimensionless Henry constant; S = transpiration rate [day"1]; F = transpiration-stream concentration factor; Dg =
gaseous diffusion coefficient [m2/d]; h is average depth of the root zone [m]; d = thickness of air boundary layer on
soil surface [m]; p = QrimRrim /QrmRrmmay be refered to as the capacity ratio; k™ = ln(2)/Xrm and kr'm =
ln(2)/Xrm are the degradation rates in the mobile and immobile phases, respectively [day"1]; Qr'm is the volumetric
immobile water content; 6rm is the volumetric mobile water content; a is a mass transfer coefficient [d"1]; Rrm and
Rr'm are the retardation factors in the mobile and immbile water regions, respectively; Drm = (Krm/6rm)^//Dg+Dz is
the effective (multiphase) dipesrion coefficient in the mobile region [m2/d]; Pr = h(v/Qrm)/Drmis the soil root-
zone Peclet number; Tr = h Rr m /(v 16 r m) is the apparent residence time in the root zone [day]; v = q + (ET 12) is
the average flux in the root zone [m/d]; and ET is evapotranspiration rate per unit area [m/d]. Equations (10-12)
describe the relationship between the leached mass emissions of organic compounds below the root zone and losses
due to biochemical degradation in the mobile and immobile phases, volatilization, and crop uptake, and the
relationship with soil and chemical properties, and related environmental factors. It extends the leached-fraction
expression of Van der Zee and Boesten (1991) to account for the effect of crop uptake, volatilization, and in
particular, nonequilibrium solute transport. Diffusive vapor losses from the immobile zone are ignored in (10).
-------
Partitioning into the vapor phase in the immobile zone, however, is allowed and can be accounted for in the
immobile-phase retardation factor Rr""= !+(/„„ pbKj+ KuimK//)/ 6rim, in which fim and Kum are the fraction of soil
bulk density and volumetric air contents in immobile water regions, respectively. And Rrm = l+(fm p\fd+ Kum KH)I
6rm, in which fm and Kum are the fraction of soil bulk density and volumetric air content in conductive pore regions
(mobile phase), respectively.
Note that that nonequilibrium transport due to diffusion into immobile-water regions has impact on the integrated
solute flux, Mr(z), only when environmental conditions are conducive for biodegradation in the immobile phase: (|>r
->• 0 as Vm ->• °° (kr'm ->• 0). In the absence of degradation in regions of stagnant water, such as intra-aggregate
pores in aggregated soils, the fraction of the solute in the immobile phase would ultimately diffuse back into the
mobile phase and leaves the root zone. The extent to which nonequilibrium transport affects the total mass available
for leaching is related to the ratio of the mass transfer rate coefficient a to the decay rate constant in the immobile
zone krim, as (11) indicates.
ETi
(a)
Precipitation
Kunn
' M°
Root
zone
q,Mr
W.T.
Aquifer
(b)
Nonequilibrium Transport
...-
/
Mpbile Soiute""
H [
\
B
Diffusive \
Transfer \
Immobile
Solute /
/
/
/
Fig. 2 (a) Illustration of a soil-aquifer system; and (b) fate and transport in mobile-immobile phase.
Eq. (10) can also be used to describe leaching and transformation in the intermediate vadose zone, assuming
negligible diffusive vapor transport and that losses are limited to biochemical degradation:
Mu(z)=Mf(h) e
r>
U
2
I T
1 i 1 u
1 1 4
i p"
ln(2)
V
V + ^u,
1
z-h
H
(13)
-------
in which 4>u is given by (11) with the associated parameters defined in the intermediate vadose zone; Pu=H(q/6um
)/Dum is the intermediate vadose soil Peclet number; and -Dum= (Kum/Qum)KHDg+Dz, the effective liquid-phase
diffusion coefficients in the intermediate vadose zone [m2/d].
The half-life parameters Xrm and Xum in (10) and (13) are understood to account for biochemical degradation. Half-
life values estimated in the field may integrate all possible loss pathways. Thus, their use in (10) and (13) may
seriously overestimate actual losses in the field. It is difficult to verify whether the values that are reported in the
literature are estimated under laboratory or field conditions. Since values for pesticides degradation half-life in
subsoils and ground water are scarce (Rao and Hornsby, 1989) and subject to uncertainty, the rational approach is to
consider "keff = lrm /[l+|a+c|>r] a random variable. In reality estimates of recharge, soil parameters, and other
environmental factors are subject to uncertainty and they display complex spatial and/or temporal variability. And a
stochastic framework, in this case, is necessary for a meaningful vulnerability assessment scheme.
UNCERTAINTY ANALYSIS
First-Order and Second-Order Analysis
Let/X) be a continuous and differentiable function of a set of N random variables X=(xi, x2, x3,.., XN), then it can be
shown by expanding/X) in a Taylor series about the mean values E(x{)= xl, E(x2)= x2, ..., E(xN)= XN , where E(*)
denotes the expectation operator, that the second-order approximation of the mean of/(X) is given by (Ang and
Tang, 1975)
f(X)=E[f(X)]=f(Xl,X2,....,xN)+-Zi:Cov(Xl,X]) ^ ' (14)
and the variance to the first-order is given by
VaV(X)} = E(f(X)-7(X)}2 = ZZ^^W!,.;) (15)
i=\i=\ dxt dxj
In which the partial and mixed-partial derivatives are evaluated at the mean vector X. In the analysis below we will
consider the following random vector X= (foc, Koc, H, X, b, 6S, \|/s, Ks). Cross correlation among these parameters is
ignored due to lack of data that support otherwise.
Uncertainty Analysis of Recharge
The mean recharge q can be approximated to the first-order by applying Eq. 14 to Eq. 9,
b
2b+3
which is (9) evaluated at the means of all random parameters. In equation (6), we have assumed yr to be
deterministic and at or greater than field capacity (-1/3 bar = 340 cm).
The first-order approximation to the varince of q is given by
(16)
dadb dH
in which partial derivatives are evaluated at X = (foc ,H,'k,b,Qs,\\is,Ks).
-------
Var(ys) and Var(b) can be obtained from the literature (Clapp and Hornberger, 1978). Var (H) is inferred from
random variations of H in time at each observation well, and an average variance is considered for each soil
hydrologic group (A,B,C,D), by averaging the variance of the wells located in each hydrologic group. The variance
of Ks however is not accounted for in this study. It should be emphasized that Var(H) is expected to vary in space
due to recharge and other effects such as boundary conditions, and that the influence of unknown pumping activities
and unknown intensive irrigation activities may result in biased estimates for Var(H).
We will compare in the application section recharge estimates based on the local equation (16) and estimates of
recharge using the SCS method (Chow et al., 1988). The latter method relies on the Hydrological groups of the soil
map, whereas the former requires the knowledge of local soil type and related properties. It is important to note that
recharge analysis on the basis of Eq. (9) does not require direct estimates of precipitation, abstractions, irrigation,
and evapotranspiration (ET). Instead, it relies on the response of the water table to the net balance in the form of
deep percolation. Its applicability therefore depends on the availability of measurements of H, where the water table
fluctuates only in response to the recharge process (i.e. precipitation and irrigation). Note that Eq. (17) relates the
variance of net recharge q directly to the temporal variability of the water table fluctuations, and not to the
variability of climatic factors, including ET. Estimation of recharge variance using the direct method would have
been difficult due to uncertainty in irrigation rates and events, and difficulties in the estimation of the variance of
evapotranspiration, due to lack of real-time measurements. Evapotranspiration (ET) is usually estimated using quasi-
empirical methods, which rely on energy balance and heat transfer concepts, and empirical crop factors. Thus,
estimation of the variance of ET requires continuous measurements in time. A modification of (9) is required for the
case of a stratified soils below the root zone and be achieved using effective soil parameters (e.g., Gelhar, 1993).
Stochastic Analysis of Pesticides Leaching
As mentioned earlier, we consider a random half-life parameter A, rather than A.e^ = V" /[l+M.+'H, which accounts for
the individual effects of all possible loss-pathway, including the parameters ^ and c|)r. This is mainly for two reasons.
First, the uncertainty in reported half-life values in the literature. And secondly, while the impact of diffusive
transfer into immobile zones can be significant, the lack of information regarding the environmental factors and
geometric properties, which are associated with or necessary for characterizing the soil mobile immobile zones,
makes it difficult to obtain meaningful estimates for the parameter 4>r . While the treatment of nonequilibrium
transport of pesticides in the soil is important, especially in aggregated soils (e.g., van Genuchten and Wierenga,
1976), it is, however, beyond the scope of this paper. We will ignore biodegradation in the intermediate-vadose soil
and limit our analysis to leaching below the biologically active root zone, at z =h. Biodegradation is likely to be
much less significant in the root zone where the likelihood of having a biotic degradation is small (Rao and
Hornsby, 1989).
Again, we make use of (14) to obtain second-order approximation of the leached mass fraction in (10):
M =M e
r ° 2 5A2 2 dr
(18)
in which the over-bar notation implies the ensemble mean and the partial derivatives are evaluated at the mean
parameters (A,, Tr) . The variance can be obtained by applying (15) to (10) evaluated at z = h,
T/ tit N ln2(2)M/
Far(Mr) = — .^ Var(Tr) + \ -=- \ Far (A,) (19)
where
l + 41n(2)—-=
-------
Var(Tr)= i Var(q) + \=pb2foc } Var(Koc) + (=pb2Koc} Var(foc) (20)
UJ ^v J U J
In Eqs. (18)-(20), we assumed that dispersion is predominantly mechanical mixing; i.e., Pr = h/aL , where aL is the
longitudinal dispersivity [m]. Thus, Pr is considered deterministic for a fixed h. Also, because the variance of ET
cannot be estimated directly for the reasons explained above, it is assumed deterministic and, thus, Var( v ) = Var (q)
follows directly from v = q + (ET 1 2) .
In the above analysis we adopt the definition of the contamination potential of a pesticide as "the likelihood of a
chemical to pass through the biologically active surface soil horizons (usually top 1 m)," (Khan and Liang, 1989).
We therefore focused our analysis above on leaching below the root zone and assumed that Mu = Mr, on the basis
that biological activity seldom extends below this zone (HUM, 1982; and Rao and Hornsby, 1991). However,
chemical reactions will further degrade some pesticieds after leaching below the root zone (Rao and Hornsby,
1991), and further analysis of (13) in the intermediate vadose soil will be required..
Pesticides concentrations in ground water are estimated using the relationship
~
g
* nB
in which Cg is the mean solute concentration in ground water [Kg/m3]; n is the aquifer porosity, and B the thickness
of contaminant plume [m]. The variance of Cg, thus, is given by
(22)
in which n and B are assumed to be constant.
APPLICATION WITH ARCVIEW CIS
Case Study: Mid- Atlantic Coastal Plain Agricultural Watershed
The study site (Locust Grove) is located in Kent County, Maryland, and covers approximately an area of 120 km2
(Fig. 3). The site is composed of two agricultural watersheds that are drained by the Morgan Creek and the
Chesterville Branch (Fig. 3). Both streams drain in to the Chester River, which lies along the southern boundary of
Kent County, Maryland, and along the northeastern shore of the Chesapeake Bay. Agriculture dominates the
landscape in the study area and most of the land there is used to grow corn and soybeans in an annual rotation with
winter wheat. Riparian forested areas are marginal and restricted to strips along the two creeks and some of the land
is used to grow ornamental trees. The soil in the study area is predominantly silt loam and to a lesser extent ranged
from loam to sandy and gravely loams. The are is underlain by the surficial (Columbia) aquifer, which consists of
sand and gravel of fluvial origin.
Arc View GIS database is developed for the soil and land use properties using the soil survey of Kent County,
Maryland (United States Department of Agriculture, Soil Conservation Service in cooperation with the Maryland
Agricultural Experiment Station and the Kent Soil Conservation District, Issued January 1982). The database is
complemented by soil hydraulic properties, including organic carbon fractions from available literature (Ravi and
Johnson, 1991, and Corselet al, 1988).
The pesticide chosen in this study is dicamba, which is among the most detected residues in the Delmarva Peninsula,
-------
and commonly used for Wheat-barley-alfalfa and Soybeans (Koterba et al., 1993). The results presented below
(GIS maps) are not based on actual agricultural land use, rather, they represent the vulnerability of ground water to
dicamba applied to soils that are suitable for growing winter wheat and small grains. A single application of 2.5
Ib/acre is assumed in the simulated vulnerability maps. An average half-life value of 16 days and a standard
deviation of 7 days are assumed for the short-lived dicamba (an Internet source provided by the ARS and the Natural
Resources Conservation Service, currently available at the following web site address:
http://www.arsusda.gov/rsml/ppdb2.html). The mean and variance of the organic carbon partition coefficient Koc
were estimated from dicamba database, using the same source: Koc =7.6 I/Kg and ^Var(Koc) =7.27 I/Kg. The
plume thickness B in ground water is assumed to be 5 m and the porosity of surficial Columbia aquifer porosity n in
the study site is considered to be 0.3 (Reilly et al., 1994).
Fig. 3 Geographical Location and detailed map of Locust Grove Study Site
Discussion of Results
Data for depths to water table are based on water levels in wells recorded in 1991, 1992 (source: Maryland
Geological Survey), 1997, and 1998 (source: US Geological Survey, Maryland). Recorded water levels in Wells KE
Be 49 and KE Be 62 (Fig. 3) from October 97 to June 98 are shown in Figs 4(a) and 4(b). In general, observed water
table levels indicated an increasing trend toward a steady level, from the month of October to the month of June,
which indicates a recharge activity. Figure 5(a) compares the theoretical recharge (Eq. 9, assuming \\is within
standard deviations reported by Clapp and Hornberger (1978)) and estimates based on soil Hydrologic Classes
using the SCS abstraction method and the method of Thornthwaite and Mather (1957). The SCS estimates are based
on monthly record of precipitation obtained from Chestertown weather station. The estimates of q using Eq. (9) were
made at several observation wells, some of which are shown on Fig. 3, based on depths to water table averaged over
the same period of November to April, from 1991, 1992, and 1998 water table records. Those values of q in
observation wells that are located in a given soil Hydrologic Class (A, B, C, and D) were averaged and compared in
Fig. 5(a) to the value estimated for the same soil Class using the SCS abstraction and Thornthwaite and Mather
-------
(1957) method (will be referred to as SCS hereafter). The apparent correspondence between the empirical and
theoretical ~q in Fig. 5(a) indicates that local recharge estimated by (9) averaged over soil Hydrologic Classes may
be consistent with the SCS recharge estimates. Thus, Eq. 9 may be used to estimate recharge at the local scale and as
a function of local soil hydraulic properties and depth to the water table. And this may hold for relatively shallow
water table in order to allow for meaningful time averages on which (9) is based.
(a)
(b)
Willl€Bj62
f
0)
o
125
12
11.5
11
1Q5
10
95
9
85
8
Ag97 Sp97 Nv97 JfenSB Feb93
Time
10
95
i 9
o
J 7
5 65
6
JLnSB JU-93
IXhA97 Jtn-93 FebSB
"fine
JLnSB
Fig.4 Records of Depth to water table in: (a) Well # KEBe59; and (b) Well #KEBe62
Figure 5(b) shows that soils in Chesterville Branch appear to be primarily class B, whereas they are primarily class
C in Morgan Creek. Thus, cumulative infiltration based on the SCS method is expected to be less in the latter due to
greater potential for surface runoff, i.e., more poorly drained. Low drainage capacity is attributed to a combination
of shallow depth to the water table and lower infiltration capacity (Philips and Bachman, 1996). This is consistent
with (9), which indicates smaller recharge, thereby greater potential for surface runoff for shallow depths to water
table, H. It also predicts less q, thus, greater surface runoff for increasingly finer soils due to smaller saturated
hydraulic conductivity Ks and lower infiltration capacity. As expected, Fig. 6(a) shows that recharge estimates using
the SCS abstractions method corresponds with the soil pattern according to the Hydrological Classification in Fig.
5(b). The estimates of mean and variance of recharge (q and Var(q)) using Eqs. (16) and (17) are shown in Figs
7(a) and 7(b). Comparison with Fig. 6(a) shows that the SCS recharge estimates are over estimated by ~q (Eq. 16),
however, both predicted a general pattern of relatively greater recharge on the east side and north of the Chesterville
Branch and around the north (headwaters) of the western branch of Morgan Creek. Also, both predicted greater
values around the north (headwaters) of the eastern branch of Morgan Creek. It is not clear if uncertainty in the
kriged estimates of depths to the water table were responsible for the relatively lower estimate of
-------
reliable strategy for the estimation of H, which is currently under investigation. In light of uncertainty in the vertical
soil variability and lack of sufficient temporal records of water table fluctuations, further assessment and evaluation
of the recharge relationship (9) is needed.
Figures 7(c-f) show first-order approximations to the mean and variance of each of dicamba mass emission to the
water table, Mu, and the resulting concentration in ground water, Cg . Biochemical degradation is ignored in the
intermediate vadose zone, and therefore Mu = Mr. This may be considered a worst case scenario, although it is based
on a single application of dicamba of 2.5 Ib/acre, and multiple applications are not uncommon. Figure 7(e) shows
that significant portions of the study area, on the east and around Morgan Creek, displayed mean groundwater
concentrations between 120 ppb and 144 ppb. These values are close to the Lifetime Health Advisory Level (HAL)
of 200 ppb (M-g/L or 10"6 Kg/m3) set for dicamba. We emphasize that the relatively greater variances in Figs. 7(d)
and 7(f) are attributed to the large variances of the residence time in the root zone, Var(Tr) (some values where of
two orders of magnitude greater than the mean), Which are the result of relatively large variances of H and Koc
parameter. This is in spite of ignoring decay losses in the intermediate vadose zone. The coefficient of variation (CV
= standard deviation/ mean) for predicted dicamba concentration in ground water, CV(Cg) = VVar(Cg)/ Cg =
(.0783)1/2/(0.000144) = 1887 (see, Figs. 7(e) and 7(f)). This implies that, although few, the estimated mean dicamba
concentrations may be subject to very high uncertainty. It is possible that first-order approximation of the variance
highly overestimated the variance. It is also likely that uncertainty in depths to the water table and uncertainty in the
above mentioned environmental factors were a contributing factor.
(a)
(b)
HYDROLOGICAL SOIL GROUPS
0.1 0.2
q (SCS) (cm/d)
0.3
RATESOF AH &• C! |
INFILTRATION
N A
Ht-ji^m* i^-i
Fig. 3. (a) Comparison between theoretical ana JsCS method ; (b) Hydrologic Soil Classification (Adapted
from Hantushetal., 1999)
Comparison between Figs. 6(b) and 7(e) clearly shows that Cg values are higly underestimated on the basis of the
SCS recharge. In fact, the spatial distribution of Cg (Fig. 6(b)) did not reflect the spatial pattern of the recharge
distribution in Fig. 6(a), rather they showed almost uniformly low concentrations (less than 1 ppb) over the two
watersheds. This is expected because of the longer residence time in the soil due to smaller SCS recharge estimates.
Also, note that recharge estimates are averaged over the soil hydrologic class in the case of the SCS method, where
as the conceptual recharge model (16) reflects the spatial variability of the soil type. And this, as well as the
distribution of the depths to the water table, may have contributed to greater predicted mass emissions to the water
table. Although the spatial distribution of the mean of dicamba concentrations in ground water (Fig. 7(e)) reflected
-------
the same spatial pattern of the mean recharge q (Eq. 16), some of the low recharge areas however showed
unexpectedly greater residual levels. This is a manifestation of the interactions among the physical and biochemical
processes and spatial variability of the environmental factors, which affect the degradation of dicamba. Recall that
the (apparent) residence time is not only a function of infiltration but also a function of the retardation factor.
Remarkably, Var(Mu) and Var(Cg) displayed a different spatial pattern as shown in Figs. 7(d) and 7(f). This
emphasizes the notion that the combined effect of spatial variability and uncertainty of the hydrologic and
environmental factors, and not exclusively convection, is responsible for uncertainty in the predicted dicamba
residuals. Although low variances of Mu and Cg were representative of the area in general, some portions, especially
those surrounding the Morgan Creek showed greater variances. The coefficient of variation, however, remains high
over the entire area.
CONCLUSIONS
A stochastically based methodology is presented for the assessment of groundwater vulnerability to nonpoint source
pesticides in agricultural watersheds. A new relationship is presented for the estimation of seasonally averaged local
recharge in terms of soil properties and seasonal average of depth to the water table. Application to a real case study
showed a correlation between recharge estimates obtained by the theoretical equation with the net infiltration
estimated on the basis of the SCS abstraction and Thornthwaite and Mather (1957) method. Physically based mass
fraction models are also presented, which describe leaching of pesticides in the root zone and the intermediate
vadose zone. The net effect of nonequilibrium transport on residual mass emissions is shown to be dependent on
ratio of mass transfer coefficient and the first-order degradation rate in the immobile phase.
A stochastic framework is developed, using first order analysis, which describes uncertainty in the recharge and
residual mass emissions of a given pesticide below the biologically active root zone. The framework was integrated
with Arc View GIS and was applied to Mid-Atlantic coastal plain agricultural watersheds for the assessment of
potential shallow groundwater contamination by the pesticide dicamba. Soil parameters and chemical and
environmental parameters (e.g., half-life, organic carbon fraction, partition coefficient and others) are subject to
uncertainty, thus, were modeled as random variables. And recharge and groundwater vulnerability to nonpoint
source of dicamba was therefore evaluated under conditions of uncertainty. The models predicted significant mean
residual levels of dicamba in ground water, which compared close to the HAL limit. The predicted mean values,
however, should be viewed with uncertainty as the relatively large variances of the mean groundwater
concentrations indicated. Results showed that the use of the SCS abstraction methods in groundwater pollution
assessment may result in the underestimation of a potential pollution problem. The predicted mean and variance of
dicamba concentrations in ground water are not exclusively driven by advective transport, rather, they reflect in their
spatial distribution the complex interactions among the physical and biochemical processes, and related
environmental factors.
The presented conceptual framework and the application to the Locust Grove study area clearly demonstrated the
order of magnitude of expected uncertainty when fixed input parameters are considered in fate and transport models.
Future attempts should account for uncertainty in NFS groundwater vulnerability assessments.
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
Ang, Alfredo H.S. and W. H. Tang, 1975. Probability Concepts in Engineering Planning and Design. John Wiley &
Sons, Inc.
Bear J. 1972. "Dynamics of Fluids in Porous Media." American Elsevier, New York.
Campbell, G. S.. 1974. "A simple method for determining unsaturated conductivity from moisture retention data."
SoilSci., 117,311-314.
-------
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., L. A. Mulkey, M. N. Lorber, and L. B. Baskin. 1985. "The Pesticide Root Zone Model (PRZM): A
Procedure for evaluating pesticide leaching threats to groundwater." Ecological Modelling. 30, 49-69.
Clapp, R. B. and G. M. Hornberger. 1978. "Empirical Equations for Some Soil Hydraulic Properties." Water
Resource. Res. 14(4), 601-604.
Corwin, D.L.; Vaughan, P.J.; League, K. 1997. "Modeling Nonpoint Source Pollutants in the Vadose Zone with
GIS." Environmental Science & Technology. 31 (8): 2157-2175.
Hantush, M. M., Z. Zhonglong, V. Murray, and M. A. Marino et al., 1999." Leaching Models for Subsurface
Pollution Assessment in Agroecosystems." ASCE-CSCE Conference on Environmental Engineering, July
25-28, Norfolk, VA.
Hantush, M., M. R. Islam, and M. A. Marno. 1998. "Models for Soil-Aquifer Vulnerability Assessment and
Pollution Prevention." Submitted to J. Hydrol.
Hillel, D. 1982. Introduction to Soil Physics. Academic Press, Inc.
Hoogeweg, C.G., and A.G. Hornsby. 1995. "Pesticide risk assessment in a watershed using GIS." Applications of
GIS to the modeling of non-point source pollutants in the vadose zone, ASA-CSSA-SSSA Bouyoucos
Conference, Mission In, Riverside, CA.
Kleveno, J. J., K. League, and R. E. Green. 1992. "Evaluation of a pesticide mobility index: impact of recharge
variation and soil profile heterogeneity." J. Contain. Hydrol., 11, 83-99.
Khan, M.A., and T. Liang. 1989. "Mapping pesticide contamination potential." Environmental Management 13:233-
242.
Koterba, M.T.; W.S.L. Banks, and RJ. Shedlock. 1993. "Pesticides in Shallow Groundwater in the Delmarva
Peninsula." J. Environ. Qual. 22: 500-518.
League K., Corwin D. L. and Ellsworth T. 1998. "The Challenge of Predicting Nonpoint Source Pollution. "
Environmental Science & Technology / News. March 1.
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.
Rao, P.S.C., A.G. Hornsby, and R.E. Jessup. 1985. "Indices for ranking the potential for pesticide contamination of
groundwater." Soil Crop Sci. Soc. Fla. Proc., 44, 1-8.
Rao, P.S.C., and A.G. Hornsby. 1991. "Behavior of Pesticides in Soils and Water." Soil Sci. Fact Sheet, SL-40
(Revised).
Ravi, V. and J. A. Johson. 1991. VLEACH A one-dimensional Finite Difference Vadose Zone Leaching Model.
USEPA R. S. Kerr Environmental Research Laboratory.
-------
Reilly, T. E., L. N. Plummer, P. J. Phillips, and E. Busenberg. 1994. " The use of simulation and multiple
environment tracers to quantify groundwater flow ina shallow aquifer." Water Resour. Res. 30(2), 421-433.
Thornthwaite, C.W. and Mather, J.R. 1957. "Instructions and Tables of computing potential evapotranspiration and
water balance." Publications in Climatology, X(3), Drexel Institute of Technology, Centerton, NJ, 1957.
Tim, U. S., JainD. and Liao H. H. 1996. "Interactive modeling of ground-water vulnerability within a geographic
information system environment." Ground Water. 34(4): 618-718.
van der Zee, S.E.A.T.M. and J.J.T.I. Boesten. 1991. "Effects of soil heterogeneity on pesticide leaching to
groundwater." Water Resour. Res., 27(12), 3051-3063.
van Genuchten, M. T. and P. J. Wierenga. 1976. "Mass Transfer Studies in Sorbing Porous Media I. Analytical
Solutions." Soil Sci. Soc. Am. Proc. 40(4). 473-480.
Wagenet R.J. and J.L. Huston, 1986. "Predicting the Fate of Nonvolatile Pesticides in the Unsaturated zone."
J. Environ. Qual. 15 (4): 315-322.
Biography
Mohamed Hantush received his B. Sc. from Kuwait University in 1985, his MSc in 1988, and his PhD in 1993, all in
Civil Engineering from the University of California at Davis. He worked as a Postdoctoral Researcher at the
University of California at Davis. He is currently a Hydrologist at the U.S. EPA National Risk Management
Research Laboratory. He is a member of ASCE, AGU, and NGWA.
Zhonglong Zhang received his B.S. from the Hebei Agricultural University in 1986 and his M.S. from the China
Institute of Water Resources and Hydroelectric Power Research (IWHR) in 1989. He has been a research engineer
in China. He also got his Ph.D in Biosystems Engineering from the Clemson University in 1998. He is working as a
programmer analyst in ManTech Environmental Research Services Corporation at Ada, Oklahoma. He is a member
of AGU and AS AE.
Miguel A Marino is professor of hydrologic science, civil & environmental engineering, and biological &
agricultural engineering at UC Davis. He has been on the faculty at UCDavis since 1972, after receiving his Ph.D.
degree in civil engineering systems from UCLA. He is the recipient of numerous awards, including the Julian Hinds
Award and Honorary Membership in the American Society of Civil Engineers.
Enamul Hoque is a senior scientist of ManTech Environmental Technology, Inc. (919 Kerr Lab Drive, Ada, Ok
74820, (580) 436-8656). He holds Masters degree in Civil Engineering from the University of Idaho. He is
responsible for managing hydrogeologic characterization studies, ground water flow and contaminant transport
modeling, and remedial designs. His research interests include nonaqueous phase liquids, transport modeling, and
development of treatment technologies for contaminated soils.
-------
(a)
(b)
•
0- 0.00004
0.00004-0.00012
| 0.00012 -0.00034
| 0.00034-0.00171
0.00171 -0.0052
0- 0.000001
0.000001 -0.000002
0.000002 -0.000015
D.OODD1S -O.DODD34
0.000034-0.00008
Fig. 6 (a) Recharge SCS (m/d); and (b) Mean groundwaterC based on SCS recharge (Kg/m )
-------
(a)
0- 0.002
| 0.002 - 0.007
I 0.007 - D.014
j 0.014- 0.026
I 0.026 • 0.041
D- 0.00018
0.00018 -0.00063
| 0.00063 -0.00132
| 0.00132 -0.00266
0.00266 -0.00488
(c)
(d)
0- 0.00017
| 0.00017 -0.00018
I 0.00018 -0.0002
0.0002 - 0.00022
'••*
[ V,
0- 0.00002
0.00002 -0.0002
| 0.0002 - 0.00284
| 0.00284 -0.01423
| 0.01423 -0.16612
(e)
(f)
•*ff'i
''.'V
t '
jf
J X
0- 0.000107
0.000107 -0.000119
0.000119 -0.00013
0.00013 -0.000144
0- 0.00009
0.00009 -0.00045
| 0.0004S -0.00126
| 0.00126 -0.00632
0.00632 -0.07383
Fig. 7 Arc View display of (a) q (m/d) (Eq. 16); (b) Var q (m4/d2); (c) Mu (Kg); (d) Var(Mu) (Kg2); (e) Cg
(Kg/m3); and (f) Var(Cg) (Kg2/m6)
-------
Precipitation.
ET
Runoff.
Soil Surface
Fig. 1 Illustration of water balance in unsaturated zone.
-------
ROOt
Zone
ET1
Precipitation
Kunoji
M°
/'?.":... '
Mu
Aquifer
Nonequilibrium Transport
..•• ••..
"
B
\
Diffusive \
Transfer \
Fig. 2 (a) Illustration of a soil-aquifer system; and (b) fate and transport in mobile-immobile phase.
-------
(a)
(b)
125
12
£11.5
,_: 11
§10.5
2 10
fas
o 9
85
8
Nw97 J=nS6
Time
JLnSB JU-93
10
95
Willl€Bj62
Jfen-93 R&93 /ig-93 JLn93
"fine
Fig.4 Records of Depth to water table in: (a) Well # KEBe59; and (b) Well #KEBe62
-------
Fig. 3 Geographical Location and detailed map of Locust Grove Study Site
-------
(a)
(a)
HYDROLOGICAL SOIL GROUPS
0.1 0.2
q (SCS) (cm/d)
0.3
Cl I OH Mi \
•s="«» «••'<•«» «
Fig.5. (a) Comparison between theoretical and SCS method ; (b) Hydrologic Soil Classification (Adapted
fromHantushetal., 1999)
------- |