Analytical modelling of pesticide transport from the soil surface to a drinking water well

Journal

Journal of Hydrology 169 (1995) 209-228

PI

Analytical modelling of pesticide transport from the soil surface to a drinking water well W.H.J. Beltmana3*, J.J.T.I. Boestena, S.E.A.T.M. van der Zeeb aDLO Winand Staring bDepartment

Centrefor Integrated Land, Soil and Water Research, P.O. Box 125,

6700 AC Wageningen. Netherlands of Soil Science and Plant Nutrition, Wageningen Agricultural 6700 EC Wageningen. Netherlands

University,

P.O. Box 8005,

Received 3 May 1994; revision accepted 17 September 1994

Abstract Pesticide transport through the unsaturated zone was modelled with an analytical solution of the convection-dispersion equation assuming steady water flow, a linear sorption isotherm and first-order transformation kinetics. Pesticide behaviour in the saturated zone was described with an analytical solution of the mass balance equation for a cylindrical flow system assuming steady flow, no dispersion, linear sorption and first-order transformation. This simplified model for the unsaturated-saturated soil system was developed to identify the processes and parameters with the greatest impact on the fraction of applied pesticide reaching a drinking water well. Leaching from the unsaturated zone was highly sensitive to the parameters describing travel time and transformation rate. Leaching increased when heterogeneity of the soil was taken into account. Pesticide arrival in the well was only moderately sensitive to the characteristic travel time and transformation rate in the aquifer. However, this sensitivity increases if zones without pesticide application were introduced around the wells (protection zones). For representative sandy soils under average Dutch rainfall conditions, processes in the unsaturated

zone had a much larger impact on pesticide arrival in the wells than processes in the saturated zone. Protection zones reduced pesticide transport to wells substantially if their half-life was much smaller than the characteristic travel time of the pesticide in the aquifer.

1. Introduction Pesticide use in horticulture and agriculture threatens groundwater, which is a major source of drinking water in many European and North American countries. At least 20 pesticides have been detected in groundwater in West European countries, * Corresponding author 0022-1694/95/$09.50 0 1995 - Elsevier Science B.V. All rights reserved SSDI

0022-1694(94)02622-X

W.H.J. Beltman et al. 1 Journal of Hydrology 169 (1995) 209-228

210

List of symbols c cl+

D Do

fom FO F a#

FU F, F,$, H K K om L Ld mx M

P Pc Per 4 r re rF’ R, RU su s t t, T T, T!l 2) V, % x z

concentration in the liquid phase (M Lm3) concentration in well water (M Le3) hydrodynamic dispersion coefficient (L2 T-‘) molecular diffusion coefficient in water (L2 T-‘) mass fraction of organic matter in soil (M M-‘) mass fraction of load transmitted by the saturated zone (M M-‘) mass fraction of load transmitted by the saturated zone; protection zone around well (M M-‘) mass fraction of dosage leached from homogeneous column (M M-l) mass fraction of dosage leached from heterogeneous field (M M-l) mass fraction of dosage arrived in well (M M-‘) thickness of the aquifer (L) slope of the linear sorption isotherm (L’ M-t) coefficient for distribution over organic matter and water (L3 M-‘) depth of the unsaturated zone (L) dispersivity (L) mean of variate X (1) area1 mass of pesticide applied (M Le2) groundwater recharge rate (L T-‘) Peclet number (1) Peclet number for heterogeneous field (1) content of pesticide sorbed (M M-‘) radius around well (L) radius of groundwater extraction zone (L) radius of protection zone around well (L) retardation factor in the aquifer (1) retardation factor in the unsaturated zone (1) relative sensitivity of the leached fraction to the ratio of travel time and half-life (1) rate of transformation of pesticide in soil (M Lm3 T-‘) time (T) travel time of solute from radius r to the well (T) characteristic travel time of water in the aquifer (T) characteristic travel time of solute in the aquifer (T) travel time of solute in the unsaturated zone (T) flow velocity of porewater in the unsaturated zone (L T-‘) flow velocity of porewater in the saturated zone in the radial direction towards well (L T-‘) flow velocity of porewater in the saturated zone in the vertical direction (L T-‘) depth below soil surface (L) depth in the saturated zone below the phreatic surface (L)

Greek letters c effective porosity of saturated zone (1) variation coefficient (1) 11 0 volume fraction of liquid in the soil (L3 Le3) x tortuosity factor for diffusion in the liquid phase in soil (1) first-order rate coefficient for transformation of pesticide (T-l) p first-order rate coefficient for transformation of pesticide in the aquifer (T-l) IAl first-order rate coefficient for transformation of pesticide in the unsaturated zone (T- -1 1 !4l dry bulk density of the soil (M L-‘ ) P

W.H.J. Tli 7, 4 Z

Beltman et al. / Journal of Hydrology

169 (1995) 209-228

211

half-life of pesticide owing to transformation in the unsaturated zone (T) half-life of pesticide owing to transformation in the saturated zone (T) frequency of pesticide application (T-l) reduction factor for fraction arriving in well owing to a protection zone (1) piezometric head (L)

with concentrations mostly between 0.01 and 100 pg 1-l (Leistra and Boesten, 1989). The maximum admissible concentration for pesticides in drinking water is 0.1 pg 1-l in the European Community (Council of the European Communities, 1980). Hence, the concentrations monitored in groundwater have exceeded this maximum up to several times. One of the aims of pesticide registration policy is to prevent contamination of drinking water wells. This policy cannot be based on monitoring studies, because transport of pesticides from the soil surface to a pumping well can take several decades. Simulation models can be valuable tools for assessing the expected pesticide arrival in drinking water wells. Wagenet and Hutson (1986) Leonard et al. (1987) and Boesten and Van der Linden (1991) developed numerical models to assess pesticide leaching from the unsaturated zone. The models describe leaching for specific crops, soil-pesticide properties and meteorological conditions. For screening purposes, simple models were developed to rank pesticides with respect to their leaching potential, based mainly on mobility and transformation parameters of the pesticide (Jury et al., 1983, 1984: Rao et al., 1985). Such ranking models do not take into account fieldscale spatial variability of soils, which may be considerable for filtration velocity (Nielsen et al., 1973; Biggar and Nielsen, 1976) organic matter content (Jokinen, 1983; Wood et al., 1987; Boekhold et al., 1991) transformation rate (Parkin and Shelton, 1992) and other soil parameters (Beckett and Webster, 1971). Both Jury and Gruber (1989) and Van der Zee and Boesten (1991) developed analytical solutions to estimate the effect of spatial variability on pesticide leaching. Van der Zee and Boesten (1991) showed that the leached fraction increased with increasing spatial variability of soil parameters. Van Ommen (1986) and Raats (1978) used exponential travel time models to describe solute transport in the saturated zone based on advective transport only. Additionally, Duffy and Lee (1992) showed that dispersion did not affect solute transport if the maximum horizontal travel distance to a small outflow interface exceeds ten times the thickness of the aquifer. Our aim was to identify processes and parameters in the coupled unsaturatedsaturated soil system that have a large impact on the fraction of the pesticide dose that reaches the pumping well. We used simplified models that are based on the models described above. 2. Theory 2.1. Homogeneous Solute

transport

unsaturated

zone

was described

by the convection-dispersion

equation,

assuming

W.H.J. Beltman et al. / Journal

212

of Hydrology 169

(1995) 209-228

first-order transformation kinetics. The equation for one-dimensional constant bulk density (p), volume fraction of liquid (0) hydrodynamic coefficient (D) and filtration velocity (v), is given by

flow, with dispersion

where q is the mass adsorbed per mass of solid phase, t is time, c is the concentration in the liquid phase, x is the depth below the soil surface and S is the transformation rate of the pesticide. The coefficient of hydrodynamic dispersion (D) in Eq. (1) is the sum of the tortuosity-corrected molecular diffusion coefficient in water, and the convective dispersion D = ADo + Ldu

(2)

where X is the tortuosity factor, Do the molecular the dispersivity. We describe pesticide sorption sorption isotherm

fom&m~

4 =

diffusion coefficient in water and Ld to the solid phase by the linear

(q)

(3)

where the mass fraction of organic matter (f,,) and the sorption coefficient for organic matter (K,,) determine the sorbed amount. The transformation rate of pesticide (5’) is described with a first-order rate equation for dissolved and adsorbed pesticide (4)

s = P(PI + Qc)

where p is the transformation rate coefficient for the liquid phase as well as solid phase. The half-life of a pesticide is given by r = ln(2)/p. Uptake of pesticides by plants is ignored; under conditions that favour small leached fractions plant uptake is usually small (Boesten, 1991). The fraction of a pesticide dose that leaches from the unsaturated zone has been derived by Jury and Gruber (1989). We used the subscript u to denote the unsaturated column zone in the equations below. The leached fraction, Fu, from a homogeneous of length L, after application of an area1 mass, M, of pesticide is F

=

u

wx =L) M(x=O)

(5) F, =exp(-kPe{

[l _“~~“1’;‘1})

where the travel time is T,

=R,L

(6)

V

and the retardation

factor is

W.H.J. Behan

&

=

et al. / Journal of Hydrology 169 (1995) 209-228

(7)

1 + f$,K,m

The Peclet number

is defined

as

Because effects of molecular diffusion sion effects, Eq. (8) may be simplified

2.2. Heterogeneous

unsaturated

213

are usually very small compared to Pe = L/Ld.

with disper-

zone

In the field, leaching parameters are spatially variable. Such a field can be visualized as an ensemble of non-interacting stream tubes, each represented by their own set of parameters. We assume steady flow, uniform properties in the vertical direction and the absence of transverse interactions between stream tubes. If X = T,/r, E p,,R,L/v is normally distributed, the average leached fraction, F,, of a heterogeneous field is (Van der Zee and Boesten, 1991) F: =exp{-kPef[(l+$)‘-l]}

(9)

in which m, is the mean (expected) value of X. The parameter Pe/is an apparent ‘field scale’ Peclet number that accounts for the variability of X. It is related to the variation coefficient, vx, by 2 vi, =pel Valocchi variation Van der because velocity. of spatial

(10)

(1985) and Jury and Sposito (1985) showed that Eq. (10) accounts for the in travel times in a semi-infinite system in the absence of transformation. Zee and Boesten (1991) pointed out that Eq. (10) can be used for X = v, the variation in travel times is determined by the variation in filtration Their approach implies that Eq. (10) can be used to account for the effect variability on leaching with X = p,,R,L/v as well.

2.3. Saturated

zone

As regards pesticide transport to a drinking water well in a phreatic aquifer, we assumed that (1) the aquifer is homogeneous isotropic and of constant thickness, (2) the lower boundary is impermeable, (3) the well fully penetrates the aquifer and (4) water flow is steady (Fig. 1). Moreover, we ignore lowering of the groundwater table owing to pumping, which is a valid assumption for distances from the well larger than 1.5-2 times the thickness of the aquifer (Bear and Verruijt, 1987, p. 23). Furthermore, we assume that pesticide sorption in the aquifer follows a linear sorption isotherm as

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

214

Drinking water well

;;

Impermeable layer Fig. 1. Stream pattern for solutes in a homogeneous calculated with Eq. (A6).

aquifer with a fully penetrating

well. Streamlines

were

in the unsaturated zone (see Eq. (3)) and that transformation is a first-order process which occurs in the solid and liquid phases (see Eq. (4)). The water flow in this system can be derived from the Laplace equation for cylindrical co-ordinates (Bear and Verruijt, 1987; p. 390): Id

rar --(

,%

ar >

+azcl’,O

(11)

az2

with r defined as the distance from the well, z defined as the depth below the phreatic level, and + as piezometric head. Darcy’s law describes the water flow in the aquifer. We assume a hydrostatic pressure head distribution with depth (Dupuit assumption), hence the radial flow velocity towards the well, II,.,is independent of the depth, z. It can easily be shown that under steady flow conditions Eq. (11) and the boundary conditions for this system are satisfied by

r2-

r,’

(124

” = 2(cH/P)r H-z ‘uz=(EHIP)

where w, is the water flow velocity in the vertical direction in the aquifer, re is the radius of the extraction area, E is the porosity, His the aquifer thickness, and P is the groundwater recharge rate. The parameter (CHIP) is the characteristic travel time of the water in the aquifer. Dispersion can be ignored when re/H > 10 (Duffy and Lee, 1992). Hence, considering only the advective transport of solutes, the mass conservation equation for the pesticide in the cylindrical system is

aq ac 1 d(rqc) PZ+EZ+Y dr

+

a(ew,c) + Pa(fC + P4)

Y$y

= 0

(13)

The derivatives to II, and U, in the velocity terms above cancel each other out, so it

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

follows

215

that (14)

We used the subscript a to denote the saturated zone. The above equations yielded the concentration in the well for a pulse dose (F,M) of pesticide which enters the top of the aquifer (see the Appendix):

c,(t) =

in which

y&w-(pa+jt-)t]

T, was the characteristic

(15)

travel time for a solute,

R,cH T, = P

defined

by

(16)

If the transformation rate coefficient pa is set at zero, Eq. (15) implies an exponential travel time distribution for a conservative solute, as derived by Van Ommen (1986). The result equals the travel time distribution of a perfectly mixed reservoir (Levenspiel, 1972). Duffy and Lee (1992) showed that this exponential travel time model also describes the transport of non-reactive solutes in hydrologically heterogeneous aquifers. For solutes subject to heterogeneous sorption and transformation, the model still applies when the integral scales of heterogeneity are small compared with the lengths of the stream tubes. We defined F, as the fraction of the applied pesticide that eventually arrives in the well. Hence, F, could be calculated from

s cc

7VZP

c,dt

= F,m;F,,M

(17)

0

Combination

of Eqs. (15) and (17) resulted

F, = lam

in

$exp[-(/ly +$)t]df WV

‘a= (Tn/TJ

k(2)

+ 1

where 7, = ln(2)/pa. 2.4. Protection

zone around well

The arrival of pesticides in a well can be reduced which pesticide use is prohibited. A protection zone radius, rp, from the well. The boundary radius, rp, related to the travel time tp of solute which enters the

by creating protection zones in may involve the area within a of the protection zone may be aquifer at this boundary by (see

216

W.H.J. Beltman et al. / Journal of Hydrology

the Appendix,

169 (1995) 209-228

Eq. (A4)): 1

tP = T, In

[ 1 - k&J2

1

(19)

Hence, travel times shorter than tp do not contribute to pesticide transport. fraction that arrives in the well, Fa,P, was calculated analogous to Eq. (18) as Fa,P = 6

The

&exp[-(p,+$-)t]dt (20)

Fas = (Ta,7,)

mc2) + 1 exp

-KLlrJ

ln(2) +

11%

It should be noted that the definition of Fas (which was based on Eq. (17)) implies that F,,P is the fraction of the hypothetical total amount of pesticide applied, assuming no protection zone. We defined a reduction factor, cp, as the quotient F,,P/FO. Combining Eqs. (18) and (20) resulted in

2Ta W)l(~.fl) cp= l-01 5 re

(21)

[

As the relative radius of the protection zone (rP/re) approaches unity, cp approaches zero and the effectiveness of the protection zone (which is (1 - cp) x 100%) is perfect. As rP/re approaches zero, and cp approaches unity, the effectiveness is zero. 2.5. Coupled unsaturated-saturated

system

Eqs. (18) and (20) yield the fraction of the mass entering the saturated zone which arrives in the well. To assess the mass fraction related to the mass input at the soil surface (the pesticide dose), we combined Eqs. (5) and (20): F, = Fu.FaP

F,+

(22a)

=exp(-kPe{

[1 +41~(‘2~TU]‘-l})

. [~ln~2)+1]-‘exp{-[$ln(2)+1]$}

Wb)

If one wishes to account for the spatial variability has to be used instead of Fu.

of the unsaturated

zone, F, (Eq. (9))

3. Results and discussion Fig.

2(a)

shows

that

pesticide

leaching

to

the

groundwater

increased

with

W.H.J.

0

10

Beltman et al. / Journal of Hydrology

20

30

40

50

L

0

217

169 (1995) 209-228

I

I

10

I(

1”’

11

20

30

40

50

T”‘B (b)

T”‘r”

Fig. 2. (a) Fraction of pesticide leaching to groundwater, F,, as a function of the ratio of travel time and half-life, T,/r,, in the unsaturated zone, as calculated with Eq. (5), for values of Pe indicated in the figure. (b) Absolute value of relative sensitivity (ls,I) of the leached fraction, F,, to the ratio of travel time and halflife, T,/r,,, in the unsaturated zone, as calculated with Eq. (26), for values of Pe indicated in the figure.

decreasing Peclet number (Pe). Hence, pesticide leaching increased with increasing dispersion. Leached fractions decreased strongly with increasing Tu/rufor Pe values in the range between five and infinity. This decrease was strongest for Pe = 00 (piston flow). If pesticide half-lives were much longer than the travel time (i.e. small Tu/7,), leached fractions were close to lOO%, irrespective of Pe. It should be noted that the vertical axis of the figure is logarithmic and ranges from unity to 10e6. Hence, Fig. 2(a) implies that the capacity of the unsaturated zone to reduce pesticide leaching may be very large. The figure shows that the leached fraction was 100% in the extreme case Because Pe = L/Ld,Pe = 0 implies either of Pe = 0, irrespective of the value of T,/r,. that the depth of the soil system, L,is zero or that the dispersivity, Ld,is infinite. The case of L = 0 is of no interest. It is remarkable that Ld = cc led to 100% leaching for high values of Tu/ru. The reason is that the analytical solution of the leached fraction (Eq. (5)) applies to a semi-infinite system. Hence, in this system the applied pesticide dose immediately spreads out between x = 0 and x = cc if Ld = cc .This means that nearly all pesticide has passed depth x = L instantly after application. The Peclet number may characterize heterogeneity of water velocity at the column scale. However, the Peclet number may also include heterogeneity of other pesticidesoil parameters at the field scale (see Eqs. (9) and (10)). We have estimated a realistic value of Peffor heterogeneous soils. Groundwater levels in Dutch non-irrigated sandy soils are mostly around 1 m below the soil surface. A dispersivity (Ld)of 0.05 m was assessed in a field experiment by Van Ommen et al. (1989) which yields Pe = 20. A Pef of 20 relates to a coefficient of variation (CV) of 0.32 for Xin Eq. (lo), if we regard the filtration velocity u as the only source of variation. To estimate Per for a heterogeneous field with Eq. (10) we need CV values for L,R, and pL, as well. Elabd and Jury (1986) and Wood et al. (1987) studied spatial variability of pesticide sorption on a scale of 1 ha and found variation coefficients of the sorption coefficient, K, of 0.310.45. Walker and Brown (1983) and Parkin and Shelton (1992) found variation

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

218

coefficients of the transformation rate coefficient between 0.07 and 0.21 and between 0.16 and 0.40, respectively. Hence, CV values of 0.30 for ~(,LL,) and n(R,) appeared reasonable values, where n(R,) < n(K) ( see Van der Zee and Boesten, 1991). We assume that ,u”, R, and ‘u are uncorrelated, and that spatial variability of L can be ignored. The variation coefficient of ,LL~R,L/w (which equals (7’,/7,)ln(2)) was approximated using ~2(&JI~)

= r12(~u) + r1*(&) + n*(L) + n*(u)

(23)

(Van der Zee and Boesten, 1991). With V(U) = 0.32, n(L) = 0, q(R,) = 0.30 and Q(,Q,,) = 0.30 we find a variation coefficient for p,,R,L/u of 0.53, which corresponds to Per = 7.1 (Eq. (10)). Based on the above-mentioned estimation, we regard Pef = 5 as a possibly realistic lower limit for non-irrigated Dutch sandy soils. Therefore, the curves in Fig. 2(a) also apply to pesticide leaching from a heterogeneous field for Pef = 0~) , 20, five and zero, respectively. These values correspond to CV values for T,,/ru of zero, 0.31, 0.63 and infinity as calculated with Eq. (10). Thus, the curves in Fig. 2(a) imply that pesticide leaching increased with increasing heterogeneity of pesticide-soil parameters. We are interested in the sensitivity of the fraction leached, FU, to the quotient T,,/r,,. Sensitivity is the rate of change in one factor with respect to the change in another factor. For this purpose we used the concept of relative sensitivity as defined by McCaen (1973) and used by Boekhold and Van der Zee (1991). The relative sensitivity (s) is defined as (24) wheref is the dependent variable and y is the independent variable. Relative sensitivity values do not depend on the magnitude off and y, and thus provide a valid means of comparing the sensitivity of the leached fraction to different values of the ratio T,/r,,. The relative sensitivity, s,, of F, to Tu/ru is given by

(25) This resulted

in

V’ul~u)W)

s, = -

U+

kWW’el(W~u)~t

(26)

The minus sign in Eq. (26) implies that s, decreases with increasing T,,/ru. The absolute value of the relative sensitivity, sU, increased with increasing T,/r,. Fig. 2(b) shows that the absolute value of this relative sensitivity increased with increasing Pe as well as with increasing T,/r,,. The leached fraction was most sensitive to Tu/r,, for Pe = co (piston flow). Ifs, is constant it follows via integraton of Eq. (25) that constant; hence a relative sensitivity of s, F,, = b(T,17JS”, where b is an arbitrary implies that the leached fraction is proportional to (Tu/~u)“u. In general, the level of the sensitivity of F, to T,/T, was large for Pe larger than five in Fig. 2(b). The curves

W.H.J.

Beltman et al. / Journal

of Hydrology169

(1995) 209-228

219

Fig. 3. Fraction

of the pesticide entering the upper groundwater that eventually arrives in the extraction of the ratio of characteristic travel time and pesticide half-life, To/r,, of the aquifer. were calculated with Eq. (20) for values of the relative radius of the protection zone, rP/re. in the figure.

well, Fos, as a function Fractions indicated

in Fig. 2(b) may also be interpreted as representing the impact of field scale heterogeneity via Pep The sensitivity s, of the leached fraction F,, from heterogeneous fields to Tu/ru decreased with increasing heterogeneity (Pef small). The leached fraction was insensitive to Tu/ru when the variation in Tu/ru was largest (Pef = 0). The calculations in this study on the leaching from the unsaturated zone were done with a simplified model (using the analytical solution of Eq. (5)). Van der Zee and Boesten (1991) analysed results of calculations with a more sophisticated numerical model which assumes Freundlich sorption, first-order kinetics for transformation and transient water flow. The transformation rate was assumed to depend on soil moisture content, temperature and depth in the soil. Van der Zee and Boesten found that the fraction leached, as calculated with the numerical model, was accurately described by Eq. (5). This indicates that the sensitivity analysis presented here probably applies to the more sophisticated models as well. Boesten (199 1) carried out a sensitivity analysis for the numerical model (for Pe = 20) and found that sensitivity to sorption and transformation parameters generally increased with decreasing leaching levels. We obtained the same result with the analytical expression, showing that this expression can be used to assess sensitivity of pesticide leaching to the input parameters. For the saturated zone, we calculated the fraction (Fa,p)of a pesticide pulse entering the aquifer that eventually reaches the well, as a function of To/r, (Fig. 3). The figure shows that the fraction Fa,pdecreased with increasing Ta/ra. Protection zones with a relative radius (Y~/Y,) of 0.2 reduced the arrival only moderately. However, for t-,/r, = 0.5, the reduction was considerable. Fig. 3 also shows that the effect of including a protection zone increased with increasing To/r,. This may be explained qualitatively as follows. If in an aquifer transformation occurs, the contribution of a stream tube to the amount of pesticide that reaches the well decreases with increasing

220

W.H.J. Beltman et al. 1 Journal of Hydrology

0

0.2

04

0.6

169 (1995) 209-228

0.8

1

r/re

Fig. 4. Reduction factor, cp. of the amount of pesticide arriving of the relative radius of the protection zone, r,,/r,, calculated travel time and half-life, Ta/ra, are indicated in the figure.

eventually in extraction wells, as a function with Eq. (21). The ratios of characteristic

travel time. The higher the transformation rate, the more important the contribution of stream tubes with short travel times will be. Protection zones exclude those stream tubes with shortest travel times. Fap is less than unity for To/r, = 0, because Fas is defined as the fraction of a diffuse input over the entire extraction area in our model. It follows that Fap is equal to the fraction of the surface area treated with the pesticide at To/r, = 0. The capacity of protection zones to reduce the amount of pesticide that reaches drinking water wells is illustrated by Fig. 4. If the pesticide is not degraded in the aquifer (T,/r, = 0) the reduction equalled the fraction of the area covered by the protection zone. Protection zones achieved only a minor decrease in the pesticide arrival at the low ratios of T,/T,. For high values of T,/T, a substantial reduction of pesticide arrival can be achieved by preventing pesticide use in a relatively small area around the well. Therefore the higher the transformation rate, the larger was the effect of introducing a protection zone. Comparison of saturated transport (Fig. 3) with unsaturated transport (Fig. 2(a)) shows that the decrease in the leached fraction with T/r found for the saturated zone is much smaller than that found for the unsaturated zone. The curve for Pe = 5 in Fig. 2(a) shows a much steeper decline than the curve for rp/re = 0 in Fig. 3. The cause of this difference can be found in the difference in the travel time distributions. Pe = 5 in the semi-infinite system for which leaching from the unsaturated zone was derived implies that the distribution of travel times is negatively skewed around the average. Travel times close to zero are very unlikely. By contrast, the exponential travel time distribution of the saturated zone has its maximum at zero travel time. The exponential travel time distribution corresponds to the finite system of the aquifer bordered by the top of the aquifer and the well tube. Consequently, a small part of the pesticide dose enters the well immediately after leaching, which is inherent to the

W.H.J.

Beltman et al. / Journal of Hydrology

169 11995) 209-228

221

(4

‘t

Fig. 5. Mass fraction of pesticide dose that eventually arrives in extraction wells, F,,, as a function of the ratios of travel time and half-life for the unsaturated zone (r,/~“) and the saturated zone (Ta/ra). The calculations were done with Eq. (22) for the relative radius of the protection zone, rP/re, set to zero. The lines are contour lines for values of fractions, F,,, indicated in the figure: (a) entire range on logarithmic scales; (b) most vulnerable range on linear scales.

description of the system. When a protection zone is instituted, the first part of the exponential distribution is cut short, so the largest fractions are left out for entering the well. Fig. 3 shows that protection zones with a relative radius rp/r, = 0.2 decreased the fraction arriving in the well only slightly compared with rp/re = 0. In the coupled unsaturated-saturated system, the fraction of the applied dose, F,,, that eventually reaches the well was calculated as a function of both Tu/ru and To/r,. This was done for rp = 0 and Per = 5, which is assumed to be a reasonable lower limit for the unsaturated zone. The result in Fig. 5(a) shows that F,,, decreased with both increasing T,,/r,, and increasing To/r,. Fig. 5(b) shows F, more detailed for high values of F,., with linear axes. Moving along a vertical line in Fig. 5(b), the decline in T,,/ru was rapid, which indicates that F,, is sensitive to T,/r,. By contrast, the contour lines are nearly parallel to the horizontal axis, indicating that F,, is not sensitive to T,/r,. Hence, the result shown in Fig. 5(b) indicates that the fraction of the pesticide dose that reaches the well was more sensitive to the T/r ratio for the unsaturated zone than to the ratio for the saturated zone,

222

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

We tried to determine realistic values for characteristic travel times and pesticide half-lives in aquifers. The thickness of Dutch phreatic aquifers varies from 10 to 250 m (VEWIN, 1986). Porosities of sandy aquifers can vary between 0.3 and 0.45. In the Netherlands, the average precipitation surplus is about 0.3 m year-‘. With these values for aquifer thickness, porosity and precipitation surplus the characteristic travel time (7’J of non-sorbing solutes ranges from about 10 to 400 years. Drainage systems such as tile-drains and ditches may collect part of the precipitation surplus and take the water out of the groundwater extraction area. In this situation, only a fraction of the precipitation surplus will reach the groundwater. This means that the characteristic travel time may be longer. Organic matter contents of aquifer material are mostly below O.l%, except for aquifers containing peaty layers. Thus, in a worstcase analysis, sorption may be ignored. Available information in the literature indicates that half-lives in aquifer material may range from less than 1 day to more than 5 years (Bromilow et al., 1986; Boesten and Van der Pas, 1993). Therefore, the range of values of T,/ra to be expected in practice is very wide (of the order of l-100 000). Groundwater is usually extracted through a partially penetrating filter, instead of the fully penetrating filter that we assumed in our model. Travel times increase if a partial filter is taken into account, because additional transport in the vertical direction is required. The travel time increases most close to the well because horizontal travel times are short there. Van der Eem (1991) showed theoretically that vertical travel times in the aquifer are negligible for the partially penetrating wells used for public water supply in the Netherlands. Thus, we assume that a partial filter will lead to approximately the same results as those shown in Fig. 3. Another model assumption was a constant thickness of the aquifer over the extraction area. In reality, the well causes a drawdown of the groundwater table in its vicinity, resulting in a greater thickness of the unsaturated zone. Consequently, travel times in the unsaturated zone increase around the well and fractions leaching to the groundwater could decrease (Fig. 2(a)). Bear and Verruijt (1987) estimated that the drawdown takes place over a distance of 1.552 times the thickness of the aquifer. For most Dutch extractions this distance is less than 0.2 times the radius of the extraction area. In Fig. 3 we saw that the introduction of a protection zone with relative radius t-,/r, < 0.2 had little impact on the fraction arriving in the wells. Therefore, ignoring the drawdown area will have little effect on the simulated arrival of the pesticide in the wells. We have derived the equations for an idealized circular extraction area. Van Ommen (1986) and Lerner and Papatolios (1993) showed that the exponential travel time approach was also valid if the shape of the extraction area is affected by, for example, streams or background regional flow. Hence, our simplified model also applies to such cases. Duffy and Lee (1992) studied the impact of random spatial variability in an aquifer-stream setting. The concentration entering the stream was not affected by spatial variation in the input concentration at the top of the aquifer, except for cases where the correlation scale of the input was larger than 10% of the maximum horizontal travel distance (our r,). This justifies our approach of using the spatially averaged input concentration, instead of distinguishing the individual fields. Duffy and Lee also found that the mean solute arrival in the stream was not sensitive to

W.H.J.

Beltman et al. / Journal of Hydrology

169 (1995) 209-228

223

random spatial variability of the hydraulic conductivity of the aquifer. This result supports the use of our analysis also for aquifers which are randomly heterogeneous in their physical properties. The model described by Eq. (22) resulted in a mass fraction of a certain pesticide application that eventually reaches the well. However, regulating authorities are interested in concentrations rather than in a fraction of a dose. The concentration in the upper groundwater is determined by the quotient of leached mass and percolated amount of water (the precipitation surplus). The model calculates a mass fraction leached; e.g. 0.1% of a realistic dose of 1 kg ha-’ yields an area1 mass in leached of 0.1 mg mP2. This leached fraction, F,, may arrive in the groundwater the course of a number of years. However, if the pesticide is applied at a certain frequency, and a steady state has developed, the concentration in the well can be estimated from

FwM4 P

c, = -

(27)

where 4 is the frequency of application. As was described above, c,, = 0.1 mg mP3 is the admissible concentration in drinking water (Council of the European Communities, 1980). As a typical example, we assume an annual application rate of 1 kg ha-’ and a recharge rate, P, of 0.3 m year-‘. It follows from Eq. (27) that F,, is 3 x lop4 at the limit value of c,. This example shows that the very small fractions arriving in the well are especially relevant in the evaluations. For example, if F,,, is 10e4 (see Fig. 5) Eq. (27) calculates an average concentration in the well of 0.03 mg mP3, which is just below the maximum admissible concentration.

4. Conclusions The unsaturated zone may have a large capacity to reduce pesticide leaching to the groundwater. Leaching from the unsaturated zone increases with increasing dispersion. The effect of dispersion on pesticide leaching from a homogeneous field is equivalent to the effect of spatial variability on pesticide leaching from a heterogeneous field. Hence, pesticide leaching from a heterogeneous field increases with increasing spatial variability of the soil. The relative sensitivity of leaching to the ratio of travel time and half-life is greatest for small dispersion, low spatial variability and low leaching levels. The saturated zone has only a small capacity to reduce pesticide transport from the top of the aquifer to the well, compared with the capacity of the unsaturated zone to reduce transport from the soil surface to the groundwater. Reduction of pesticide arrival in drinking water wells by protection zones becomes increasingly effective for pesticide residues with increasing transformation rate in the aquifer. Considering the coupled unsaturated-saturated transport, we conclude that in most situations the fraction of the dose leaching from the unsaturated zone to the groundwater will be smaller than the fraction of the load passing the aquifer to the drinking water wells. The European Community guideline for drinking water can be

224

W.H.J. Beltman et al. / Journal of Hydrology

169 (1995) 209-228

met only when a very small fraction (e.g. 0.0001) of the pesticide dose reaches the drinking water well. To reduce pesticide arrival in drinking water wells, the main measure should be to reduce pesticide leaching, because the reduction capacity of the saturated zone is limited. If pesticides leach significantly to the groundwater, drinking water wells are vulnerable to pesticides with long half-lives in the aquifer. Under these conditions, drinking water wells are difficult to protect with protection zones. Therefore, further emphasis should be put on research assessing the transformation of pesticides in aquifers.

Acknowledgements We would like to thank J.L.W. Gielen of the Department of Mathematics, Agricultural University, Wageningen, for discussions regarding the Appendix, and M. Leistra (SC-DLO) for useful comments on the manuscript. This work was part of the Netherlands Integrated Soil Research Programme (W.B, J.B.) and partly funded by the Directorate-General Science, Research and Development of the Commission of the European Communities, via the European Community Environmental Research Programme STEP-CT90003 1 (S.v.d. Z.).

Appendix: Mathematical solution transport in saturated zone Eq. (14) is a first-order partial differential equation. The partial differential equation can be replaced by an equivalent system of ordinary differential equations (Zachmanoglou and Thoe, 1986). Those ordinary differential equations can be solved. We used the variable T defined as T = CHIP. The characteristic of Eq. (13) is given by (t(s), r(s), z(s), c(s));

d*

a=

Ra

=s

dr_3

r2 - r,2

H-z

pdz=ds

=+

- -R,,LL~c

;dc

+

+

H-z

=+

t = R,s+a Tln(r$r2)

= -Rap0

-Tln(H-z) *

230

c=o

=s+y

c = 6exp(-R,p,s)

in which a (T), p (T), y (T) and 6(M Lp3) are integration and boundary conditions t=O

=s+,B (Al)

T

ds=--Fds-

-dr=ds

2Tr

dz

+

2Tr

j

ds-

de

dt = R,ds

constants.

Using the initial

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

z=o

t>O

225

c = co

we determined the characteristic going through assumed that the pesticide front of concentration water table at t = 0. We followed the front entering the given point Y coincide with s = 0, giving

(t, r, z, c) = (0, ro, 0, co) = Y. We co enters the groundwater at the the groundwater at r = ro. We let

a=0 ,L3= Tln(r,2 - ri) y = -Tin(H)

642)

s = co Hence,

the above characteristic

through

Y is as follows:

t = R,s

(A34

Tln

6-b)

Tln

(A3c)

c = co exp(-Rapas)

(A34

Eliminating s results in a full description of the solute behaviour in the radial system. The travel time t from radius ro to the well follows from Eqs. (A3a) and (A3b): t = R,Tln

1 [ 1 - (r0/rA2

1

Eqs. (A3a) and (A3c) yield the following

(A4) relationships

between

t and z:

Wb) Hence, solute travel time to depth z does not depend on radius r (see, e.g. Raats, 198 1, Fig. 9). The velocity dz/dt at z = 0 equals (H/R,T). The travel time to the lower boundary of the aquifer is infinite. Solute velocities V, = dr/dt and V, = dz/dt can be derived from Eqs. (A4) and (A5). They are related to the velocities in Eq. (12) via V, = q/R, and V, = Q/R,. The path of a streamline starting at radius r. was found by combination of Eqs. (A3b) and (A3c):

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

226 2

2

z,ffro-’ ri

-

646)

r2

This equation can be used to draw the streamlines in the aquifer (see Fig. 1). The relation between depth and concentration was derived from Eqs. (A3c) and (A3d), giving (A7) It should be noted that c does not depend on the distance from the well, r. For a penetrating front, c(z) can only be calculated for z < zP, where zP is the penetration depth of the front to be calculated with Eq. (A5a). For z > zP, it follows that c(z) = 0. The water flow rate at the well is horizontal and does not depend on z. Therefore, the concentration in the well, c,, was simply obtained by averaging over the thickness of the aquifer: 1 =P c,, = c(z)dz Hos Substitution

of Eq. (A7) gave (A9)

The concentration in the well as a function of time was derived function of time via Eq. (A5) for zP in the above equation:

by using

z as a

cdt) (AlO) -=p,~~+l{l-exp[-(iLn+~)IlJ co The concentration in the well following a pulse input (Dirac delta function, S(t)) with mass M (M LP2) is determined by the derivative of the above equation to t:

where Eq. (17) was used to eliminate

co.

References Bear, J. and Verruijt, A., 1987. Modeling groundwater flow and pollution. Reidel, Dordrecht, 414 pp. Beckett, P.H.T. and Webster, R., 1971. Soil variability: a review. Soils Fert., 34: l-15. Biggar, J.W. and Nielsen, D.R., 1976. Spatial variability of leaching characteristics of a field soil. Water Resour. Res., 12: 78-84. Boekhold, A.E. and van der Zee, S.E.A.T.M., 1991. Long term effect of soil heterogeneity on cadmium behaviour in soil. J. Contam. Hydrol., 7: 371-390. Boekhold, A.E., van der Zee, S.E.A.T.M. and de Haan, F.A.M., 1991. Spatial patterns of cadmium contents related to soil heterogeneity. Water Air Soil Pollut., 57-58: 479-488.

W.H.J. Beltman et al. / Journal of Hydrology 169 (1995) 209-228

227

Boesten, J.J.T.I., 1991. Sensitivity analysis of a mathematical model for pesticide leaching to groundwater. Pestic. Sci., 31: 375-388. Boesten, J.J.T.I. and van der Linden, A.M.A., 1991. Modeling the influence of sorption and transformation on pesticide leaching and persistence. J. Environ. Qual., 20: 425-435. Boesten, J.J.T.I. and van der Pas, L.J.T., 1993. Transformation rate of atrazine in water-saturated sandy subsoils. In: T. Eggers (Editor), Quantitative Approaches in Weed and Herbicide Research and their Practical Application, Proc. 8th EWRS Symposium, Braunschweig, 1993. EWRS, Wageningen, pp. 381-387. Bromilow, R.H., Briggs, G.G., Williams, M.R., Smelt, J.H., Tuinstra, L.G.M.Th. and Traag, W.A., 1986. The role of ferrous ions in the rapid degradation of oxamyl, methomyl and aldicarb in anaerobic soils. Pestic. Sci., 17: 535-547. Council of the European Communities, 1980. Directive of the Council of 15 July 1980 relating to the quality of water intended for human consumption. Off. J. Eur. Commun., L229(30 Aug.): 1 l-29. Duffy, C.J. and Lee, D.-H., 1992. Base flow response from nonpoint source contamination: simulated spatial variability in source, structure, and initial condition. Water Resour. Res., 28: 9055914. Elabd, H. and Jury, W.A., 1986. Spatial variability of pesticide adsorption parameters. Environ. Sci. Technol., 20: 256-260. Jokinen, R., 1983. Variability of topsoil properties at the southern coast of Finland and the number of soil samples needed for the estimation of soil properties. J. Sci. Agric. Sot. Finl., 55: 109-l 17. Jury, W.A. and Gruber, J. 1989. A stochastic analysis of the influence of soil and climatic variability on the estimate of pesticide groundwater pollution potential. Water Resour. Res, 25: 246552474. Jury, W.A. and Sposito, G., 1985. Field calibration and validation of solute transport models for the unsaturated zone. Soil Sci. Sot. Am. J., 49: 1331-1341. Jury, W.A., Spencer, W.F. and Farmer, W.J., 1983. Behaviour assessment model for trace organics in soil. I, Model description. J. Environ. Qual., 12: 558-563. Jury, W.A., Spencer, W.F. and Farmer, W.J., 1984. Behaviour assessment model for trace organics in soil, II, Chemical classification parameter sensitivity. J. Environ. Qual., 13: 567-572. Leistra, M. and Boesten, J.J.T.I., 1989. Pesticide contamination of groundwater in Western Europe. Agric. Ecosyst. Environ., 26: 369-389. Leonard, R.A., Knisel, W.G. and Still, D.A., 1987. GLEAMS: groundwater loading effects of agricultural management systems. Trans. ASAE, 30: 1403-1418. Lerner, D.N. and Papatolios, K.T., 1993. A simple analytical approach for predicting nitrate concentrations in pumped ground water. Ground Water, 31: 370-375. Levenspiel, O., 1972. Chemical Reaction Engineering, 2nd edn. Wiley. New York, 566 pp. McCaen, 1973. The role of sensitivity analysis in hydrologic modeling. J. Hydrol., 18: 37-53. Nielsen, D.R., Biggar, J.W. and Erh, K.T., 1973. Spatial variability of field measured soil-water properties. Hilgardia, 42: 2155221. Parkin, T.B. and Shelton, D.R., 1992. Spatial and temporal variability of carbofuran degradation in soil. J. Environ. Qual., 21: 6722678. Raats, P.A.C., 1978. Convective transport of solutes by steady flows II. Specific flow problems. Agric. Water Manage,, 1: 219-232. Raats, P.A.C., 1981. Residence times of water and solute within and below the root zone. Agric. Water Manage., 4: 63-82. Rao, P.S.C., Hornsby, A.G. and Jessup, R.E., 1985. Indices for ranking the potential for pesticide contamination of groundwater. Proc. Soil Crop Sci. Sot. Fla., 44: 1-8. Valocchi, A.J., 1985. Validity of the local equilibrium assumption for modeling sorbing solute transport through homogeneous soils. Water Resour. Res., 21: 8088820. Van der Eem, J.P., 1991. Pesticides and phreatic groundwater extractions. KIWA-Report SW0 90.312, KIWA, Nieuwegein, Netherlands, 39 pp, (in Dutch). Van der Zee, S.E.A.T.M. and Boesten, J.J.T.I., 1991. Elects of soil heterogeneity on pesticide leaching to groundwater. Water Resour. Res., 27: 3051-3063. Van Ommen, H.C., 1986. Influence of diffuse sources of contamination on the quality of outflowing groundwater including non-equilibrium adsorption and decomposition. J. Hydrol., 88: 79-95.

228

W.H.J. Beltman et al. 1 Journal of Hydrology

169 (1995) 209-228

Van Ommen, H.C., van Genuchten, M.Th., van der Molen, W.H., Dijksma, R. and Hulshof, J., 1989. Experimental and theoretical analysis of solute transport from a diffuse source of pollution. J. Hydrol., 105: 225-251. VEWIN, 1986. Provincial review of resources and sources. VEWIN, Rijswijk, Netherlands (in Dutch). Wagenet, R.J. and Hutson, J.L., 1986. Predicting the fate of nonvolatile pesticides in the unsaturated zone. J. Environ. Qual., 15: 315-322. Walker, A. and Brown, P., 1983. Spatial variability in herbicide degradation rates and residues in soil. Crop. Prot., 2: 17-25. Wood, L.S., Scott, H.D., Marx, D.B. and Lavy, L.T., 1987. Variability in sorption coefficients of metolachlor on a Captina silt loam. J. Environ. Qual., 16: 251-256. Zachmanoglou, E.C. and Thoe, D.W., 1986. Introduction to Partial Differential Equations with Applications. Dover Publ., New York, 391 pp.