Semidiscrete pesticide transport modeling and application

Journal of Hydrology 285 (2004) 19–40 www.elsevier.com/locate/jhydrol

Semidiscrete pesticide transport modeling and application Xuefeng Chua,*, Miguel A. Marin˜ob a

Annis Water Resources Institute, Grand Valley State University, 740 West Shoreline Drive, Muskegon, MI 49441, USA b Department of Land, Air, and Water Resources and Department of Civil and Environmental Engineering, University of California, Davis, CA 95616, USA Received 21 January 2003; accepted 25 July 2003

Abstract A time-continuous and space-discrete method is proposed for three-phase (dissolved, adsorbed, and vapor phases) pesticide transport modeling in the vadose zone and a generalized semidiscrete solution is derived under conditions of heterogeneous media, unsteady flow fields, and space-time-dependent physical and biochemical processes concerning pesticide environmental fate. The developed model also takes into account pesticide runoff and erosion and pesticide transport in the plant canopy zone. The model is able to deal with various pesticide application methods commonly used in practice, including over-canopy, undercanopy, or any combined foliar and soil surface spray, as well as soil-incorporated applications. The hybrid semidiscrete solution method, formulated herein, incorporates analytical and numerical methodologies into a uniform, flexible modeling framework, which makes the model suitable for either screening-level investigations as a lumped model or detailed studies as a distributed model. Furthermore, the semidiscrete pesticide transport model is applied in the Orestimba Creek Basin, California, for evaluating the vulnerability of the hydrosystem to diazinon contamination. In accordance with the observations, the simulation demonstrates that diazinon exposure levels in the creek frequently exceed criteria for aquatic life, and many peak pulses even exceed the water-quality standard for human health. In the subsurface environment, however, high concentrations of diazinon are generally limited within the shallow soil. It is also found that the magnitude and timing of pesticide application and rainfall/irrigation dominate exposure levels of diazinon in both subsurface and surface environments. q 2003 Elsevier B.V. All rights reserved. Keywords: Semidiscrete method; Transport modeling; Pesticide contamination; Vadose zone; Surface water

1. Introduction Due to frequent detection of pesticides and their breakdown products in groundwater and surface water, much attention has been paid to pesticideinduced non-point source (NPS) contamination of * Corresponding author. Tel.: þ1-616-331-3987; fax: þ 1-616331-3864. E-mail address: [email protected] (X. Chu). 0022-1694/$ - see front matter q 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2003.07.004

surface and subsurface environments. Undoubtedly, enhanced ability to simulate pesticide transport and transformation and to predict the exposure levels of pesticide residues in ground and surface water would significantly benefit protection of the vulnerable hydrosystem and circumvention of further deterioration in water quality. Considerable efforts have been made to simulate pesticide fate and transport in the subsurface environment and a number of models, including both analytical and

20

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

numerical models, have been developed during the last three decades. Van Genuchten et al. (1974) described an analytical model for predicting pesticide movement in porous media. Jury et al. (1983, 1987) developed analytical models for screening pesticide fate and transport in soils under steady-state soil water flow conditions. Jury and Gruber (1989) developed an analytical model on the basis of a probability density function of the residual mass fraction of pesticide remaining after transport through the biologically active zone in order to specifically identify the relative role of soil and climatic variability on pesticide leaching. Beltman et al. (1995, 1996) derived an analytical solution for a coupled unsaturated –saturated system to determine the pesticide contamination in wells with a continuous input (1995) and certain application frequency (1996). Hantush and Marin˜o (1996) proposed a physically based analytical model for long-term prediction of pesticide contamination in a coupled soil-aquifer system. Complete mixing was assumed for the crop root zone and intermediate vadose zone while in the aquifer a two-dimensional advection– dispersion model was developed. Their model was also tested by Chu et al. (2000) by comparing against a numerical counterpart that linked PRZM2 and MT3D. The comparison indicated the limitations of the assumption of complete mixing. Furthermore, Hantush et al. (2000) extended their analytical modeling work to non-equilibrium solute transport in a dual-porous medium to analyze the effect of local dispersion in mobile– immobile phases on the pesticide leaching potential. More recently, Hantush et al. (2002) developed mass fraction models for multi-phase transport and fate of agricultural pollutants in tworegion structured or aggregated soils under steady unsaturated flow conditions. In analytical pesticide transport models (such as the aforementioned models), steady water flow and time-invariant parameters are commonly assumed. To derive closed-form solutions, boundary conditions (e.g. pesticide loading) are often simplified. In practice, however, pesticides can be applied in many ways (over-canopy or soil surface spray, and soil-incorporated application). Also, regardless of the method of application, pesticides are potentially subject to surface runoff and erosion. Existing analytical

models generally fail to account for these important factors and simulate the related processes. Thus, analytical pesticide transport models are only suitable for screening purposes. Since the early 1980s, significant advancement in numerical modeling has been achieved. The representative numerical models for pesticide transport simulation in the subsurface environment include PRZM (Pesticide Root Zone Model, Carsel et al., 1984; Mullins et al., 1993; Carsel et al., 1998), RUSTIC (Risk of Unsaturated/Saturated Transport and Transformation of Chemical Concentrations, Dean et al., 1989), CREAMS (Chemicals, Runoff, and Erosion from Agricultural Management Systems, Knisel, 1980), GLEAMS (Groundwater Loading Effects of Agricultural Management Systems, Leonard et al., 1987; Knisel et al., 1993), LEACHM (Leaching Estimation and Chemistry Model, Wagenet and Hutson, 1989), AGNPS (Agricultural Non-Point Source Pollution Model, Young et al., 1994), and HSPF (Hydrological Simulation Program—FORTRAN, Bicknell et al., 1996). These models have a variety of complexities and many of them have undergone a series of upgrades and improvements. Several of them have been widely used for simulating fate and transport of agricultural pesticides in soils and overland flow under different conditions (Carsel et al., 1986; Jones et al., 1986; Melancon et al., 1986; Pennell et al., 1990; Parrish et al., 1992; Varshney et al., 1993; Spurlock, 1998). Generally, numerical models are capable of dealing with heterogeneous media, transient flow, and dynamic biochemical kinetics. Yet, intensive data requirements, computation load, and numerical dispersion and oscillations are common concerns. Alternatively, a semidiscrete approach can be useful for solving pesticide transport problems, in which the space domain is discretized while keeping the time domain continuous. Consequently, the original second-order partial differential equation (PDE) that governs pesticide transport is converted into a system of first-order ordinary differential equations (ODEs) by using a numerical method (finite differences or finite elements). Then, the ODE system is analytically solved to yield the final closed-form solution. The semidiscrete approach has been applied to contaminant transport modeling in porous media

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

(Guymon, 1970; Guymon et al., 1970; Nalluswami et al., 1972; Willis, 1979; Hwang et al., 1984; Umari and Gorelick, 1986; Kamra et al., 1991a,b; Skaggs and Barry, 1996). This solution method overcomes the limitations to the analytical solution engendered by heterogeneous media and irregular boundary geometry (Hwang et al., 1984). Since the solution is continuous and exact in time, no time-marching is necessary for temporal computations, which significantly minimizes numerical error accumulations. Thus, solute concentrations at any spatial point can be calculated without passing through all intermediate time steps (Hwang et al., 1984; Umari and Gorelick, 1986). This technique eliminates artificial smearing and oscillations in the numerical solution induced by time discretization and the relevant schemes, such as implicit or explicit schemes. Compared to the standard numerical approach, the semidiscrete method often yields smaller truncation errors for long-term predictions (Kamra et al., 1991a,b). However, a steady-state water flow field and time-invariant parameters are commonly assumed in the aforementioned studies. To the best of our knowledge, this solution method has never been applied to three-phase pesticide transport in variably saturated porous media with time- and space-variant parameters that represent heterogeneous media, unsteady flow fields, and complex biochemical processes. In this study, a hybrid time-continuous and space-discrete (TC-SD) semidiscrete method is proposed for pesticide transport modeling, which takes advantage of both analytical and numerical methods. The overall objective of this study is to develop a pesticide transport model using the TCSD semidiscrete solution method that is able to deal with both time- and space-varying physical and biochemical processes related to one-dimensional (1D) vertical water flow and three-phase pesticide transport in the vadose zone. The model also simulates pesticide runoff and erosion and accounts for various pesticide application methods. In particular, over-canopy spray of pesticides is taken into account by analyzing pesticide fate in the canopy zone. Due to use of the semidiscrete solution method, the resulting model possesses a flexible structure that enables one to use it as either a lumped model for screening purposes or a distributed model for detailed studies. To deal with

21

a regional scale problem, the model has been programmed for section-based simulations. That is, a watershed is divided into small sections and 1D transport modeling is performed for each section. Furthermore, the model is tested by applying it to the Orestimba Creek Basin, California, for evaluating diazinon environmental fate in subsurface and surface water systems.

2. Description of the methods 2.1. Pesticide transport processes in the plant canopy – vadose zone system Following the release of pesticides into the environment, they can be volatized into the atmosphere, washed off the target fields into adjacent rivers/streams due to runoff and erosion, and leached down through the vadose zone into underlying aquifers. Pesticides exist in three phases (dissolved, adsorbed, and vapor phases) and undergo a set of physical and biochemical processes affected by numerous factors such as pesticide properties, climatic conditions, hydrogeological settings, agricultural practices, etc. This study focuses on simulating pesticide transport in the plant canopy –vadose zone system. The conceptual model structure is schematically shown in Fig. 1. In addition to the plant canopy zone, the entire vadose zone is divided into the surface zone, plant root zone, and deep vadose zone so as to account for distinct physical and biochemical features

Fig. 1. Pesticide transport in the plant canopy–vadose zone system.

22

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40



 ›C › ›ðKH CÞ › ›C ¼ aDg uR uD l þ ›z ›z ›t ›z ›z › ðqCÞ 2 ðrr þ re þ ru þ rd þ ru ÞC 2 ›z þ Mðz; tÞ ð2Þ in which R ¼ 1 þ ðrKd þ aKH Þ=u

ð3Þ

rd ¼ uRks

ð4Þ

ru ¼

Fig. 2. Pesticide processes in the vadose zone.

along the unsaturated profile. Major processes modeled for these three zones are depicted in Fig. 2. They include advection, dispersion, degradation, sorption, partitioning between dissolved and vapor phases, plant root uptake, volatilization, runoff, and erosion. A 1D physically based compartmental model is used for simulating water flow along the soil profile. Surface runoff and erosion are modeled using the soil conservation service curve number (SCS-CN) method and modified universal soil loss equation (MUSLE), respectively (Appendix B). 2.2. Mathematical expressions for pesticide transport modeling By assuming linear equilibrium sorption, linear equilibrium liquid –vapor partitioning, and first-order decay, the second-order PDE governing 1D threephase pesticide transport in the vadose zone can be expressed as

› ½ðu þ rKd þ aKH ÞC ›t

 › › › ›C aDg ðKH CÞ þ uDl ¼ ›z ›z ›z ›z › ðqCÞ 2 ðrr þ re þ ru ÞC 2 ›z 2 ðu þ rKd þ aKH Þks C þ Mðz; tÞ

ð1Þ

›ðuRÞ ›t

ð5Þ

where C is the concentration of the dissolved-phase pesticide [M/L3]; u; the volumetric water content [L3/L3]; a; the volumetric air content [L3/L3]; r; the bulk density [M/L3]; Kd ; the distribution coefficient [L3/M]; KH ; the dimensionless Henry’s law constant; Dg ; the diffusion coefficient of the vaporphase pesticide [L2/T]; Dl ; the dispersion coefficient of the dissolved-phase pesticide [L2/T]; q; the water flux [L/T]; ks ; the first-order degradation rate [1/T]; rr ; the pesticide runoff rate [1/T]; re ; the pesticide erosion rate [1/T]; ru ; the pesticide root uptake rate [1/T]; M; the pesticide loading term [M/L3/T]; R is the retardation factor. Note that all parameters, such as u; a; q; ks ; Kd ; KH ; Dg ; Dl ; R; rr ; re ; and ru ; can be functions of spatial coordinate z and time t: Determination of the pesticide loading term M depends on application methods, such as overcanopy spray, soil surface spray, and soil-incorporated application. Particularly, pesticide fate in the canopy zone is analyzed in this study so as to take into account over-canopy spray of pesticides. The relevant details are shown in Appendix A. Expressions for pesticide runoff rate, rr ; erosion rate, re ; plant root uptake rate, ru are given by Eqs. (B5), (B9) and (B14), respectively, in Appendix B. Note that Eq. (1) is a general expression. For different zones, however, it may involve different terms. For example, the governing equation for the plant root zone does not include pesticide runoff and erosion, which are special processes in the thin surface zone. The diffusion coefficient for vapor-phase pesticide transport in porous media can be expressed as (Jury et al., 1991; Simunek

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

et al., 1998) Dg ¼ jg ðaÞDag ¼ ða7=3 =n2 ÞDag

ð6Þ

where jg ðaÞ is the gas tortuosity factor; n; porosity [L3/L3], and Dag is the binary diffusion coefficient of the vapor-phase contaminant in free air [L2/T]. Similarly, the hydrodynamic dispersion coefficient (including molecular diffusion and mechanical dispersion) for dissolved-phase pesticide transport in porous media is given by (Bear, 1972; Simunek et al., 1998) lql Dl ¼ Dp þ aL lvl ¼ jl ðuÞDw l þ aL u lql ð7Þ u p where D is the molecular diffusion coefficient of the dissolved-phase pesticide in porous media [L2/T]; jl ðuÞ; the liquid tortuosity factor; Dwl ; the binary diffusion coefficient of the dissolved-phase pesticide in water [L2/T]; aL ; the longitudinal dispersivity [L]; v is the pore water velocity [L/T]. Thus, the final governing equation can be written as
 ›Cðz;tÞ › g ›ðKH ðz;tÞCðz;tÞÞ ¼ E ðz;tÞ uðz;tÞRðz;tÞ ›z ›t ›z
 › l ›Cðz;tÞ › E ðz;tÞ þ 2 ½qðz;tÞCðz;tÞ ›z ›z ›z 2rðz;tÞCðz;tÞþMðz;tÞ ð8Þ ¼ ðu7=3 =n2 ÞDw l þ aL

Eg ðz; tÞ ¼ aDg ¼ ða10=3 =n2 ÞDag l

E ðz; tÞ ¼ uDl ¼ ðu

10=3

=n

2

ÞDw l

ð9Þ þ aL lql

rðz; tÞ ¼ rr þ re þ ru þ rd þ ru

ð10Þ ð11Þ

2.3. Initial and boundary conditions 2.3.1. Initial condition Cðz; t0 Þ ¼ C0 ðzÞ

the upper and lower boundaries on the pesticide transport along a profile of variably saturated porous media. Note that volatilization of the vapor-phase contaminant through the top boundary can be approximated by an air-boundary layer model (Jury et al., 1983; 1991, refer to Appendix B.3 for details). Thus, the top boundary can be expressed as For t . 0 and z ¼ 0 : 2Eg

ð12Þ

where t0 denotes the initial time and C0 ðzÞ is the corresponding initial concentration at location z: 2.3.2. Boundary conditions Generally, Cauchy-type boundary conditions (BCs) can be used to describe the effects of

›ðKH CÞ ›C 2 El þ qC ›z ›z

¼ qCin 2

Da ðKH C 2 Cga Þ d

ð13Þ

By assuming Cga ¼ 0 and ›ðKH CÞ=›z ¼ KH ›C=›z at the top boundary, Eq. (13) can be simplified as 2ET

›C þ ðq þ Pv ÞC ¼ f0 ›z

ð14Þ

in which Pv is given by Eq. (B13) in Appendix B, and ET and f0 are expressed as ET ¼ Eg KH þ El f0 ¼ qCin þ

Da a C ¼ qCin d g

ð15Þ ð16Þ

Similarly, the bottom boundary can be written as for t . 0 and z ¼ L : 2ET

in which

23

›C þ qC ¼ qCout ›z

ð17Þ

3. Semidiscrete solution method To solve the second-order PDE (Eq. (8)) together with the initial and boundary conditions (Eqs. (12) – (17)) using a hybrid semidiscrete approach, the space domain, including the surface zone, plant root zone and deep vadose zone, is discretized into N cells. By applying finite differences or finite elements, the original PDE can be transformed into a system of ODE as shown in the following matrix form _ ¼ AðtÞCðtÞ þ MðtÞ CðtÞ

ð18Þ

where AðtÞ is the coefficient matrix; CðtÞ; the _ concentration vector; CðtÞ; the derivative vector of the concentration; MðtÞ is the pesticide loading vector. As an example, if a two-point upwind finite difference scheme is applied, the coefficient AðtÞ is a tridiagonal

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

24

matrix and the transformed ODE system can be explicitly expressed as: 3 3 2 2 a1;1 a1;2 C_ 1 ðtÞ 7 7 6 6 7 7 6 6 7 6 C_ 2 ðtÞ 7 6 a2;1 a2;2 a2;3 7 7 6 6 7 7 6 6 7 7 6 6 7 6 C_ 3 ðtÞ 7 6 a a a 3;2 3;3 3;4 7 7 ¼6 6 7 7 6 6 7 7 6 6 7 6 .. 7 6 .. 7 6 . 7 6 . 7 7 6 6 5 5 4 4 aN;N21 aN;N C_ N ðtÞ 3 2 3 2 C1 ðtÞ m1 ðtÞ 7 6 7 6 7 6 7 6 6 C2 ðtÞ 7 6 m2 ðtÞ 7 7 6 7 6 7 6 7 6 6 C ðtÞ 7 6 m ðtÞ 7 7þ6 3 7 3 6 ð19Þ 7 6 7 6 7 6 7 6 6 . 7 6 . 7 6 . 7 6 . 7 6 . 7 6 . 7 5 4 5 4 CN ðtÞ

mN ðtÞ

Specifically, the semidiscrete expression for cell i can be written as C_ i ðtÞ ¼ ai;i21 ðtÞCi21 ðtÞ þ ai;i ðtÞCi ðtÞ þ ai;iþ1 ðtÞCiþ1 ðtÞ þ mi ðtÞ in which g l ðK ðEg þ Ei21 Þ þ Eil þ Ei21 Þ ai;i21 ¼ H;i21 i ui Ri Li ðLi þ Li21 Þ 2qi21 ui Ri ðLi þ Li21 Þ g l ðKH;i ðEig þ Eiþ1 Þ þ Eil þ Eiþ1 Þ ai;i ¼ 2 ui Ri Li ðLi þ Liþ1 Þ þ

ð20Þ

2qi r 2 i ui Ri ðLi þ Li21 Þ ui Ri

g l ðKH;iþ1 ðEig þ Eiþ1 Þ þ Eil þ Eiþ1 Þ ai;iþ1 ¼ ui Ri Li ðLi þ Liþ1 Þ

mi ¼

Mi u i Ri

a1;1 ¼ 2

KH;1 ðE1g þ E2g Þ þ E1l þ E2l u1 R1 L1 ðL1 þ L2 Þ

! KH;1 E1g þ E1l q1 r 2 a 2 1 þ u 1 R1 L 1 0 u 1 R1 u1 R1 L21 ! KH;1 E1g þ E1l M1 q1 fp þ þ m1 ¼ u1 R 1 u 1 R1 L 1 0 u1 R1 L21

ð21Þ

ð25Þ

ð26Þ

in which

a0 ¼ f0p ¼

2L1 ðq1 þ Pv Þ 2ðKH;1 E1g þ E1l Þ þ L1 ðq1 þ Pv Þ

ð27Þ

2L1 f0 þ E1l Þ þ L1 ðq1 þ Pv Þ

ð28Þ

2ðKH;1 E1g

Similarly, applying the bottom BC Eq. (17), the coefficient entry for the lower boundary node, aN;N should be changed as

aN;N ¼ 2 2

g l KH;N ðENg þ EN21 Þ þ ENl þ EN21 uN RN LN ðLN þ LN21 Þ

2qN r 2 N uN RN ðLN þ LN21 Þ u N RN

ð29Þ

The solution to the non-homogeneous time-variant state equation system (18) is given by ðt CðtÞ ¼ Fðt; t0 ÞCðt0 Þ þ Fðt; tÞMðtÞdt ð30Þ t0

g l ðK ðEg þ Ei21 Þ þ Eil þ Ei21 Þ 2 H;i i ui Ri Li ðLi þ Li21 Þ

2

the coefficient entries, a1;1 and m1 ; should be changed as:

ð22Þ ð23Þ

ð24Þ

where Li is the discretized space increment of cell i: For the upper boundary node ði ¼ 1Þ; after taking into account the upper boundary condition (Eq. (13)),

in which the state transition matrix F is expressed as (Friedly, 1972) ( j21 ðt1 X ðt Fðt; t0 Þ ¼ lim I þ Aðt1 Þ Aðt2 Þ· · · j!1



ð tk t0

k¼0

t0

t0

Aðtkþ1 Þdtkþ1 · · ·dt1



ðt ðt Aðt1 Þdt1 þ Aðt1 Þ ¼Iþ t t0
ð t1 0  Aðt2 Þdt2 dt1 þ · · · 

ð31Þ

t0

where I is a unit matrix (i.e. ui;i ¼ 1 for all i; ui;j ¼ 0 for all i and j where i – j).

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

Eq. (30) represents the closed-form solution of the ODE system (i.e. TC-SD solution to the original PDE). If the input data are time-discrete values (e.g. daily data), the non-homogeneous ODE system (Eq. (18)) can be simplified as a set of piecewise homogeneous problems as follows (note that the coefficient matrix A is independent of time within each time interval):

For continuous pesticide loading MðtÞ ¼ m

ð40Þ

_ ¼ ACðtÞ þ m CðtÞ

ð41Þ

The solution is CðtÞ ¼ eAðt2t0;i Þ Cðt0;i Þ þ Since

_ ¼ ACðtÞ þ MðtÞ; CðtÞ ð32Þ t [ ½t0;i ; t0;i þ Ti ; i ¼ 1; …; NT in which t0;i ¼ t0 þ

iX 21

ð33Þ

Tj

where t0;i is the initial time for interval i; Ti ; the lasting time for interval i; NT, the total number of time intervals; t0 is the initial simulation time. Also, since parameters u and R are constant in each time interval, we have

›ðuRÞ ¼0 ð34Þ ›t Thus, the state transition matrix and the analytical solution to Eq. (32) in time interval ½t0;i ; t0;i þ Ti  can be, respectively, simplified as A2 ðt 2 t0;i Þ2 þ ··· Fðt; t0;i Þ ¼ I þ Aðt 2 t0;i Þ þ 2! ru ¼

¼

eAt ¼ I þ At þ

1 X Ak ðt 2 t0;i Þk ¼ eAðt2t0;i Þ k! k¼0

CðtÞ ¼ eAðt2t0;i Þ Cðt0;i Þ þ

ðt

eAðt2tÞ MðtÞdt

ð35Þ ð36Þ

t0;i

Further, expressions for the semidiscrete equation (32) and the corresponding solution (Eq. (36)) can be given for the following two distinct loading cases. For instantaneous pesticide loading (Dirac input) MðtÞ ¼ mdðt 2 t0;i Þ ð37Þ _ ¼ ACðtÞ þ mdðt 2 t0;i Þ CðtÞ

ðt

eAðt2tÞ m dt

ð42Þ

t0;i

1 X A2 t2 A3 t3 Ak t k þ þ ··· ¼ 2! 3! k! k¼0

ð43Þ

# 1 ðt X Ak ðt 2 tÞk dt m e m dt ¼ k! t0;i t0;i k¼0 " # 1 X Ak ðt ¼ ðt 2 tÞk dt m k! t 0;i k¼0 " # 1 k X A ðt 2 t0;i Þkþ1 m ¼ ðk þ 1Þ! k¼0 " # 1 X Ak ðt 2 t0;i Þkþ1 21 ¼A A m ðk þ 1Þ! k¼0 " # 1 X Akþ1 ðt 2 t0;i Þkþ1 21 ¼A m ðk þ 1Þ! k¼0 " A2 ðt 2 t0;i Þ2 21 ¼A Aðt 2 t0;i Þ þ 2! # 3 3 A ðt 2 t0;i Þ þ þ · ·· m 3!

ðt

j¼1

25

Aðt2tÞ

"

¼ A21 ½eAðt2t0;i Þ 2 Im

ð44Þ

Thus, the final concentration vector is given by j k CðtÞ ¼ eAðt2t0;i Þ Cðt0;i Þ þ A21 eAðt2t0;i Þ 2 I m

ð45Þ

4. Introduction to a case study in the Orestimba Creek Basin, California

ð38Þ 4.1. Background of pesticide contamination

The solution is given by CðtÞ ¼ eAðt2t0;i Þ ½Cðt0;i Þ þ m

ð39Þ

where m is the pesticide loading vector and dðtÞ is the Dirac-delta function.

Agriculture benefits substantially from use of a variety of pesticides. However, pesticide-induced contamination of the vulnerable hydrosystem has become a critical environmental issue. According to

26

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

8200-sample analyses conducted by the National Water Quality Assessment (NAWQA) program during 1992– 1996, 58 pesticides were detected at least once at or above 0.01 mg/l in both ground and surface water in 20 of the nation’s major hydrologic basins (USGS, 1998). Over the past two decades, pesticides or their transformation products have been detected in ground waters of more than 43 states. At least 143 pesticides and 21 transformation products have been detected (Barbash, 1995). Also, low levels of pesticides have been widespread in the nation’s surface waters for several decades. In many streams, some pesticides frequently exceed water quality criteria seasonally in each year although their annual average concentrations seldom exceed regulatory standards for drinking water. Three insecticides including diazinon, carbaryl, and chlorpyrifos were most frequently detected in 58 nationwide rivers and streams (Larson et al., 1999). In California, pesticides and their breakdown products have been found in 3845 wells, and 75 pesticides or their breakdown

products have been detected in 44 counties throughout California (Pease et al., 1995). For the surface water, 49 pesticides were detected in the San Joaquin River and three subbasins. The concentrations of seven pesticides exceeded the criteria for aquatic life and 40% of those exceedances were attributed solely to diazinon (Dubrovsky et al., 1998). Diazinon, one of the two primary dormant season organophosphate (OP) insecticides, is generally applied with an oil mixture in orchards in the San Joaquin Valley. Yet coincidence of such dormant-season applications with seasonal rainfall events provides a high potential for pesticides to wash off target fields into nearby surface water bodies due to runoff and erosion. Dormant spray pesticides enter the estuary in distinct pulses at concentrations that are toxic to invertebrates (Kuivila and Foe, 1995). The California Department of Pesticide Regulation has been conducting monitoring studies on acute and chronic toxicity at the Orestimba Creek site and in the San Joaquin River at Vernalis since the early 1990s. Diazinon was found at high

Fig. 3. Regional map of California’s Central Valley and location of the study area.

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

27

levels in the surface water, which is toxic to test organisms at both sites. In addition, diazinon applied during spring and summer in the Orestimba Creek Basin was also washed off the target fields and eventually discharged into the San Joaquin River with irrigation return flow (Domagalski, 1997). Thus, the Orestimba Creek Basin was selected in this case study and the semidiscrete pesticide transport model was applied for evaluating diazinon contamination in subsurface and surface environments. In particular, section-based modeling was programmed, which enables one to deal with this basin-scale problem even though only 1D simulation is performed for each small section. 4.2. Description of the study area As shown in Fig. 3, the study area is located in the Central Valley of California, which covers the western parts of Stanislaus and Merced counties and comprises three small subbasins in the San Joaquin Valley: Orestimba Creek Basin (ORE), Central California Irrigation District (CCID), and Coast Ranges (CR). The ORE, CCID, and CR subbasins include 13 sections (33,670,000 m2), 22 sections (56,980,000 m2), and 164 sections (424,760,000 m2), respectively (note that the area of a section is 1 mile2 or 2,590,000 m2). The study area in the valley floor is largely agricultural land. Diazinon is applied for agricultural pest control in both crop growth periods and dormant seasons. The temporal variations of the diazinon applications over the entire area (daily values) and the spatial distribution of the total amount of diazinon applied in each section during 1996– 1997 are shown in Figs. 4a and b, respectively (based on data provided by the California Department of Pesticide Regulation). The diazinon-applied crop fields comprise 6 sections in the ORE, 15 sections in the CCID, and 3 sections in the CR. Clearly, the heaviest diazinon applications in 1996– 1997 are located in the CCID subbasin and the intensive diazinon applications in crop growth periods and dormant seasons occurred in June and January, respectively. Almond and walnut are two primary orchard crops, and tomatoes, beans, and alfalfa are major vegetable, field, or pasture crops in the 24 diazinon-applied sections. Over-canopy spray is the major diazinon

Fig. 4. Diazinon applications in the study area in 1996–1997.

application method for this orchard-dominated study area. The crop growth periods generally range from April to September. Irrigation is essential for all crops and under-canopy furrow irrigation is commonly adopted during most of the growth periods of the crops. Except for the mountainous and foothill areas in the western part of the Coast Ranges, the study domain in the valley floor consists of confluent alluvial fans with higher clay content and hence lower permeability. According to the soil survey report by Ferrari and McElhiney (1997), clay loam, loam, and sandy loam are three major soil types in the study area. The depth of sediments increases from west to east. The water table is generally deeper than 10 m in the crop-planted areas. The precipitation primarily falls in November through March and rainfall data at Newman and Gilroy Stations in 1996– 1997 are used for the simulation. Orestimba Creek is an ephemeral stream. The water flow is dominated by rainfall-induced surface runoff from the three subbasins in the winter and it also receives irrigation return water in the spring and summer

28

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

(Panshin et al., 1998; Domagalski, 1997). As a tributary, Orestimba Creek eventually discharges into the San Joaquin River (Fig. 3).

5. Diazinon transport modeling in the Orestimba Creek Basin 5.1. Basic information and model parameters

Table 1 Soil hydraulic parameters Soil

ur

us

ns

Ks (m/day)

r (g/cm3)

FC

Sandy loam Loam Clay loam References

0.065 0.078 0.095 –a

0.41 0.43 0.41 –a

1.89 1.56 1.31 –a

1.061 0.2496 0.0624 –a

1.49 1.45 1.40 –b

0.20 0.22 0.34 –b

a b

As mentioned previously, the study area includes a total of 199 sections, 24 of which are diazinon-applied sections. The simulation period ranges from January 1, 1996 to December 31, 1997 (731 days). The semidiscrete pesticide transport model is applied to evaluate the diazinon environmental fate in the basin. The section-based modeling structure particularly facilitates basin-wide simulations. Meanwhile, spatial variations in soil types, crops, climatic conditions, pesticide applications, etc. over the basin can be taken into account in the model. In this basin-scale application study, the surface runoff and erosion are calculated over the entire study area (199 sections). However, the pesticide transport modeling in the vadose zone is conducted only for the diazinonapplied sections (24 sections). In the model, pesticide transport in the overland flow and along stream channels is not physically simulated. Instead, lumped regression models, to be discussed next, are used to describe contributions of both water and pesticide residues from each section to the outlet of the basin. Due to the low permeability of soils, the transport simulation indicates that diazinon residues primarily concentrate within the shallow soil (less than 1 m) and hence the maximum simulation depth of the vadose zone for all diazinon-applied sections is 2.52 m (the thickness of the thin surface zone is 0.02 m and the discretized space increment for the remaining 2.5-m soil profile is 0.05 m). The plant canopy zone intercepts rainfall/over-canopy irrigation water and the plant interception storage capacity is 0.08 –0.3 cm (Singh, 1992; Mullins et al., 1993). The runoff curve numbers (CN) under average antecedent moisture condition (AMC-II) are 75 – 80 for crop growth periods and 80 – 85 for no-crop periods. In the model, runoff curve numbers are automatically calculated based on simulated soil moisture conditions. Table 1 shows soil hydraulic parameters,

Simunek et al. (1998). Mullins et al. (1993).

including residual water content, ur ; saturated water content, us ; soil hydraulic characteristics parameter, ns ; saturated hydraulic conductivity, Ks ; bulk density, r; field capacity, FC for different soils. Table 2 lists the major parameters concerning diazinon properties. All model parameters are prepared for each diazinon-applied section. Many of them (e.g. soil water content, flow velocity, diazinon decay rate, distribution coefficient) are also varying in time and along the layered soil profiles. Detailed spatial and temporal distributions of the soil water content and flow velocity for each section are provided by the soil water flow model and ranges of some variable parameters (e.g. decay rate and distribution coefficient) are given in Table 2. Different values of these parameters are used for the surface zone, plant root zone, and deep vadose zone. Within each zone, the same parameters are used for all cells. It is assumed in the simulation that the initial diazinon concentrations in the vadose zone are zero for all sections. Implementation of the section-based modeling not only provides spatial and temporal distributions of diazinon concentration in the vadose zones of the simulated sections, but also the quantity of surface runoff from all sections and the corresponding diazinon concentrations, which will be further used to estimate diazinon contaminations in Orestimba Creek. 5.2. Computation of stream discharges and diazinon concentrations at the outlet of Orestimba Creek To quantify pesticide exposure levels in the surface water and evaluate the vulnerability of Orestimba Creek to diazinon contamination, the relationship of the precipitation and stream flow is analyzed and multiple linear regression models are used to compute

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

29

Table 2 Major parameters on diazinon properties Parameter

Value

Refs

Distribution coefficient, Kd Henry’s constant, KH First-order decay rate in soil, ks First-order decay/volatilization rate in the canopy zone, kc Diffusion coefficient in free air, Dag Diffusion coefficient in water, Dw l Longitudinal dispersivity, aL Water solubility, Sw Air boundary layer thickness, d Organic carbon content, foc Log of octanol–water partitioning coefficient, log Kow

0.65–3.23 cm3/g 5.0 £ 1025 (dimensionless) 0.01–0.02 day21 0.10–0.30 day21 0.43 m2/day 0.000043 m2/day 0.05 m 40.0 mg/l 0.005 m 0.001– 0.005 g/g 3.02 (cm3/g for Kow )

– a,b –a –a – a,c –d –d – b,e –f –d –b –a

a b c d e f

Mullins et al., 1993. Computed/assumed. Knisel, 1980. Jury et al., 1991. Spitz and Moreno, 1996. Kidd and James, 1991.

the total water discharge, pesticide quantity and concentrations at the outlet of Orestimba Creek (at River Road). Finally, we are able to evaluate the loading potential of diazinon from the Orestimba Creek Basin to the San Joaquin River. A relationship between observed discharges at the outlet of the basin and rainfall events can be determined by using a multiple linear regression model Ri ¼ R0 þ b0 Pi þ b1 Pi21 þ · · · þ bK Pi2K þ 1i ;

ð46Þ

i ¼ 1; 2; …; Nt If the intercept R0 is zero, Eq. (46) can be rewritten as X K Ri ¼ bk Pi2k þ 1i ; i ¼ 1; 2; …; Nt ð47Þ

the observed response Ri from the fitted value R^ i (least squares method). This process is performed by utilizing the Microsoft IMSL Stat Library and the computation is conducted separately for different climatic zones (ORE/CCID and CR). The resulting coefficients reflect simple relationships between stresses (precipitation) and responses of the system (flow at the outlet). Due to lack of the detailed information on actual irrigation management and practices, the same sets of regression coefficients are used for both rainfall and irrigation events in this study. By using the above regression relationships, the water and pesticide contributions of diazinon-applied section j to the outlet at time i can be written as

k¼0

where Ri is the observed discharge at the outlet of the basin at time i (dependent variables or the responses to the rain events); Pi2k ; the rainfall at time i 2 k (independent variables); bk ; the regression coefficients; R0 ; the base flow of the discharge at the outlet (intercept); 1i ; the independently distributed normal error with mean zero and variance s2 for time i; Nt ; the total time period number; K is the number of observations. The regression coefficients are then determined by minimizing the sum of squares of the deviations of

K X

Rji ¼

ajk Qji2k ;

k¼0

ð48Þ

i ¼ 1; 2; …; Nt ; j ¼ 1; 2; …; Nsp

Wij ¼

K X

j ajk Qji2k Ci2k ;

k¼0

i ¼ 1; 2; …; Nt ; j ¼ 1; 2; …; Nsp

ð49Þ

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

30

Table 3 Weighting factors for the ORE/CCID and CR subbasins k

0

1

2

3

4

5

6

Total

ak (ORE/CCID) ak (CR)

0.1924 0.2856

0.2406 0.3305

0.1662 0.1317

0.1073 0.0562

0.0937 0.0436

0.0735 0.1453

0.1264 0.0071

1.0000 1.0000

in which

Nsp K X X

b ak ¼ K k ; X bk

k ¼ 0; 1; 2; …; K

ð50Þ

Ciout

W ¼ i ¼ Ri

k¼0

where Qji2k is the simulated runoff for diazinonapplied section j at time i 2 k; Rji ; the runoff water contribution from diazinon-applied section j to the discharge at the outlet at time i; Wij ; the pesticide runoff contribution from section j to the diazinon discharge at the outlet at time i; Nsp is the number of the total diazinon-applied sections; ajk is the weighting factor of contribution of section j at time k: The computed weighting factors for the two climatic zones are shown in Table 3. In addition to the diazinon-applied sections, the other sections without diazinon application also contribute their surface runoff to the stream discharges at the outlet of the basin. Their water contribution can be expressed as Rji ¼

K X

ajk Qji2k ;

ð51Þ

k¼0

i ¼ 1; 2; …; Nt ; j ¼ 1; 2; …; Nnsp where Nnsp is the number of sections without diazinon application. Thus, the total amount of water Ri ; total mass of dissolved-phase pesticide Wi ; and the corresponding overall concentration Ciout at the outlet of Orestimba Creek at time i can be estimated, respectively, by using the relations: Nsp þNnsp

Ri ¼

Wi ¼

X

K X

j¼1

k¼0

Nsp K X X j¼1 k¼0

ajk Qji2k ;

j ajk Qji2k Ci2k ;

i ¼ 1; 2; …; Nt

i ¼ 1; 2; …; Nt

ð52Þ

ð53Þ

j ajk Qji2k Ci2k

j¼1 k¼0 Nsp þNnsp

X

K X

j¼1

k¼0

;

ajk Qji2k

ð54Þ

i ¼ 1; 2; …; Nt 6. Analysis of results 6.1. Diazinon contamination in soils Temporal distributions of the simulated diazinon concentration in the soil depth of 0.01 – 0.145 m (shallow soil) and 0.395 – 0.545 m (deeper soil) for two typical sections (M08S08E02 and M08S08E15, Fig. 3) are shown in Figs. 5 and 6. In the shallow soil (Figs. 5b and 6c), occurrence of diazinon peak concentrations depends on the timing of diazinon application and rainfall/irrigation. Generally, the peaks of diazinon concentration in the top soil coincide with dormant season applications in the winter due to frequent, intensive rainfall events. For instance, diazinon applications in section M08S08E15 in January 1997 resulted in the highest peak concentration. In crop-growth periods (spring and summer), the model simulation indicates that diazinon peaks generally correspond to the application events. However, the peak concentrations may not be in direct proportion to the applied diazinon amount. Clearly, if diazinon is applied by over-canopy spray and irrigation water is applied under-canopy, a large portion of diazinon may never reach the soil surface due to lack of carrying water and high degradation/ volatilization rate in the plant canopy. Thus, it can be concluded from the simulation results that pesticide application, irrigation, and irrigation methods (overcanopy or under-canopy) have significant influence on

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

31

Fig. 5. Diazinon concentrations in soil for section M08S08E02.

pesticide levels in the shallow soil during crop-growth periods. Furthermore, it can be inferred that overcanopy irrigation methods present a much higher risk of pesticide contamination to both subsurface and surface environments during crop-growth periods than under-canopy methods. Figs. 6 also indicates that although the amount of diazinon applied in June 1996 is much higher than that in January 1997 (Fig. 6b), the highest peak of diazinon in the shallow soil is located in January (Fig. 6c) due to the heavy

rainfall (Fig. 6a). In contrast, a combination of overcanopy pesticide spray and under-canopy irrigation method is commonly used in this orchard-dominated section during summer months, which tends to reduce the diazinon exposure levels in the shallow soil. In the deeper soil, however, a different pattern is observed in the model concerning temporal distributions of diazinon concentration as shown in Figs. 5c and 6d. Unlike in the shallow soil, major diazinon peaks uniformly occur at the end of January or early

32

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

Fig. 6. Diazinon concentrations in soil for section M08S08E15.

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

33

February, just following the heavy rainfall season although the simulated sections have distinct diazinon application schedules (Figs. 5a and 6b). The direct influence of diazinon applications on the simulated pesticide distributions in the deeper soil almost disappears. This occurs because vast infiltrating water and deep percolation are primarily induced by the heavy rainfall. Thus, the model suggests that diazinon concentrations in the deeper soil are mainly controlled by deep percolation generated by heavy rainfall in the winter. The simulation results also indicate that with increase in the depth, diazinon concentrations in soil decrease rapidly and diazinon exposure levels are negligible for the soil deeper than 1 m. No obvious threat of diazinon contamination to the underlying groundwater can be inferred for the study area. It is evident from Fig. 4 that heavy diazinon applications were in June of 1996 and 1997 and dormant-season sprays mainly occurred in January 1997. Thus, three times (June 30, 1996; April 30, 1997; and June 30, 1997) are selected herein to examine potential effects of the intensive applications on the subsurface environment by showing the distribution of diazinon concentrations in the shallow soil over the study area. Figs. 7a – c display the spatial variability of simulated diazinon concentrations at the depth of 0.01 m at the three selected times, respectively. The simulated concentrations of diazinon on April 30, 1997 are the highest (Fig. 7b). This can be attributed to dormant-season applications and heavy rainfall, which washes off pesticide residues and moves them into the soil. Due to a significantly large amount of diazinon that was applied in June, 1996 (Fig. 4), the concentrations on June 30, 1996 are also very high. Moreover, it can be observed that simulated diazinon concentrations in the shallow soil in the CCID subbasin are generally higher than those in the other areas, which is consistent with the spatial distribution of the diazinon applications. 6.2. Diazinon contamination in Orestimba Creek

Fig. 7. Spatial variability of diazinon concentrations over the study area.

Figs. 8a and b show the comparison of the simulated and observed water discharge and diazinon concentrations, respectively, at the outlet of Orestimba Creek (River Road) (note that the observed concentrations only range from May 23, 1996 to April

30, 1997). Generally, diazinon occurrence in Orestimba Creek is determined by the availability of both surface runoff water and contaminant diazinon. Fig. 8b indicates that there are two major diazinon peaks at the end of July 1996 and the end

34

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

Fig. 9. Overall simulated water and pesticide contributions of the three subbasins at the outlet of Orestimba Creek in 1996– 1997.

Fig. 8. Comparison of the simulated and observed discharge and diazinon concentrations of Orestimba Creek at River Road.

of April 1997, respectively, and both peaks are in dry seasons, when the stream discharge is very low (Fig. 8a). From Fig. 4, it can be concluded that vast diazinon applied in the CCID in June 1996 and a certain amount of irrigation return flow led to the first peak at the outlet of the Creek in July 1996. The irrigation return flow potentially carries a certain amount of pesticide residues into nearby surface water bodies in the primary diazinon-applied summer months (e.g. June and July). Because the water flow of Orestimba Creek is very low during the summer months, concentrations of diazinon in the Creek can be very high even for a relatively small amount of pesticide loading. Hence, irrigation can play an important role during seasons without precipitation although the irrigation return flow is limited. The second peak of the diazinon concentration in April 1997 is attributed to the dormantseason and spring applications in 1997. The irrigation return flow in the relevant sections also plays an important role in the transport of diazinon from crop fields to the stream channel and eventually to the outlet of the basin. The simulated diazinon concentrations in Orestimba Creek appear to be over-predicted for some crop growing months (June – August, 1996). This can probably be

attributed to the fact that we did not take into account the details and complexity of actual irrigation practices in the modeling. A prototype study conducted by the USGS showed high daily variability of diazinon concentrations during irrigation periods, which was driven by differences in irrigation frequency and water management (Domagalski, 1997). In the winter, dormant-season applications did not induce high concentrations of diazinon in the Creek in 1996 –1997. It can be attributed to two factors. First, a huge amount of the surface runoff from the broad diazinon-free Coast Ranges significantly dilutes the diazinon concentrations in the Creek. Diazinon concentrations tend to decrease rapidly as the flow from the CR subbasin increases. This phenomenon has been also observed in the USGS monitoring work (Domagalski et al., 1997). As shown in Fig. 9, the overall water contribution of the CR subbasin to the discharge at the outlet of Orestimba Creek in 1996 – 1997 is as high as 94% while the diazinon mass contribution is only 6%. This is an important reason for lower exposure levels of diazinon in the Creek in the winter of 1996 and 1997. Next, the amount of diazinon applied during the dormant seasons and the timing of the applications and rainfall events also affect the concentrations in the stream. Compared to the applications in crop-growth periods, less diazinon was applied during the dormant seasons of 1996 – 1997 (Fig. 4). 6.3. Water quality and management strategies Fig. 8b shows that diazinon levels in Orestimba Creek frequently exceed the water-quality standard

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

for aquatic life (8 £ 1025 mg/l, Larson et al., 1999). Worst of all, diazinon concentrations even exceed the human-health criterion (6 £ 1024 mg/l, Larson et al., 1999) and such high peak pulses may last about half a month. In the subsurface environment, however, high concentrations of diazinon are limited to the shallow soil (, 1 m). No obvious threat to the underlying groundwater system can be inferred from the simulation results due to the low permeability of the soils and the large thickness of the unsaturated zone (. 10 m). This finding resulted from the modeling is consistent with the conclusion obtained from USGS monitoring studies that organophosphate pesticides have not been detected in groundwater in the western San Joaquin Valley (Domagalski et al., 1997). The pesticide transport modeling and analyses of the simulated results, conducted herein, indicate that integrated water – pesticide –crop management can be an effective way to minimize pesticide exposure levels in surface and subsurface environments. Specifically, pesticiderelated water quality control can be achieved by adopting the following strategies: (1) During winter rainfall seasons, reducing the rainfall-induced runoff by various crop field management and adjusting the timing of dormant season applications and rainfall (try to increase their time intervals); (2) in the spring and summer dry seasons (crop-growth periods), selecting appropriate irrigation methods and scheduling that potentially reduce both surface runoff (irrigation return flow) and deep percolation while meeting the water need for crop growth; and (3) avoiding the over-canopy irrigation following over-canopy sprays of OP pesticides. In addition, modeling results indicate that as high as 66% of diazinon residues in Orestimba Creek come from the CCID subbasin (Fig. 9). Thus, the CCID is the major area for water quality control.

7. Conclusions A TC-SD hybrid solution method was proposed and a mathematical model was developed for simulating three-phase pesticide transport in the vadose zone. The semidiscrete model accounted for both spatial and temporal variability in all parameters that represent heterogeneous porous media,

35

unsteady flow, as well as space-dependent and timevariant physical and biochemical processes associated with pesticide transport and transformation. The primary physical and biochemical processes simulated in the model include advection, diffusion/dispersion, sorption, volatilization, partitioning between dissolved and vapor phases, degradation, plant root uptake, runoff, erosion, etc. occurring in the surface zone, plant root zone, and deep vadose zone. Pesticide fate in the plant canopy zone was also analyzed. The model is capable of handling various pesticide application methods commonly used in practice (e.g. over-canopy, under-canopy spray, and soil-incorporated applications). The TC-SD semidiscrete method is essentially a two-step solution approach. The first step involves spatial discretization and numerical approximation while the second step focuses on solving the resultant first-order ODE system. For the first step, different numerical schemes of varying order accuracy (first-order or high-order) and mathematical features (linear or non-linear) can be used. In addition, decoupling the spatial and temporal issues into two steps makes it possible for us to use highorder interpolation in the discretized space domain and then to employ high-order ODE solvers to solve the resulting system of ODE (Le Veque, 1992). The flexible model structure, resulting from the use of the hybrid semidiscrete solution method, particularly improves the applicability of the model. It consequently can be used as a lumped model (like a multi-zone model or box model) for screening purposes if data are sparse and quick decisions are needed. In this case, the model provides exact analytical solutions in time for each zone/box (zoneaveraged pesticide concentrations). On the other hand, the model can also be extended to a complex distributed model to account for detailed spatial and temporal variability in porous media, water flow, and pesticide transport and environmental fate if enough data are available. In the case study, the semidiscrete pesticide transport model was applied to the Orestimba Creek Basin to evaluate the environmental fate and transport of diazinon in subsurface and surface water systems. The performance of the model was tested using real-world data and the simulated concentrations of diazinon in the Creek

36

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

were compared with the measured values. However, due to lack of data, a similar comparison between measured and simulated concentrations of diazinon along soil profiles was not conducted. It was demonstrated that the magnitude and combined timing of pesticide application and rainfall dominate exposure levels of diazinon residues in both subsurface and surface environments during the winter rainy season. In crop-growth periods (spring and summer), occurrence of high concentrations of diazinon in surface water bodies and shallow soils primarily depends on diazinon application and irrigation (amount, timing, and methods). The model suggests that diazinon-induced contamination in the subsurface environment is only limited to the top soil (, 1 m) and no obvious threat to the groundwater can be inferred in the study area. In the surface water, however, diazinon exposure levels may exceed the water-quality standards for both aquatic life and human health. These modeling conclusions regarding diazinon distributions and contamination potential are consistent with the USGS findings resulted from long-term monitoring studies. Finally, integrated water – pesticide – crop management is suggested as a strategy to minimize pesticide contamination in both surface and subsurface environments.

Acknowledgements We would like to thank Graham Fogg and Gerald Orlob for their insightful review of the original manuscript. We also thank Frank Spurlock and Craig Nordmark (CalEPA, DPR) and Joseph Domagalski and Peter Dileanis (USGS) for providing valuable information and field data that have been used for model testing in this study. This work was supported by the US Environmental Protection Agency (USEPA) (Grant No. R819658) Center for Ecological Health Research at the University of California, Davis, the UC Toxic Substances Research and Teaching Program, and the UC Pesticide Management Research and Extension Project. Although the information in this document has been funded in part by the USEPA, it may not necessarily reflect the views of the Agency and no official endorsement should be inferred.

Appendix A. Pesticide application and fate in the plant canopy zone Pesticides can be applied in various methods, which lead to distinct environmental pathways and different spatial and temporal distributions of pesticide residues in the hydrosystem. For over-canopy spray application method, it is necessary to analyze the pesticide fate in the plant canopy zone. Generally, over-canopy pesticide spray is partitioned into two parts: canopy application ðMc Þ and direct soil application ðMds Þ (ignoring the portion of pesticide that is directly volatized into the atmosphere during applications) Mc;t ¼ st Mt

ðA1Þ

Mds;t ¼ ð1 2 st ÞMt

ðA2Þ

where Mt is the total pesticide application in time period t [M/L2] and st is the partitioning ratio of pesticide application between the plant canopy and soil surface (0 , st , 1 for combined foliar and soil surface over-canopy pesticide spray, st ¼ 1 for fully over-canopy pesticide spray, and st ¼ 0 for undercanopy soil surface application or soil-incorporated application). The pesticide application partitioning ratio st varies with a number of factors, such as pesticide application methods, actual pesticide application procedures, crop types, crop growth seasons, wind, etc. In this manner, the model is capable of dealing with any combined pesticide applications between the plant canopy and soil. The canopy-applied pesticide ðMc Þ can be further washed off the plant canopy zone to the soil surface (referred to as canopy-soil application Mcs ) due to the throughfall (referred to as canopy-soil rainfall/irrigation Pcs =Ics ) and the pesticide in the canopy zone is also subject to decay/volatilization. The overall mass balance of pesticide in the plant canopy zone can be expressed as dMSTG ðtÞ ¼ mc ðtÞ 2 mcs ðtÞ 2 kc ðtÞMSTG ðtÞ dt

ðA3Þ

where MSTG ðtÞ is the pesticide mass storage in the plant canopy zone [M/L2]; mc ðtÞ; the canopy pesticide application rate [M/L 2/T]; mcs ðtÞ; the canopysoil pesticide application rate [M/L2/T]; kc ðtÞ is

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

the first-order pesticide decay/volatilization rate in the plant canopy zone [1/T]. For time increment Dt; the pesticide mass storage at time t is given by

however, it is assumed that the pesticide is evenly distributed within the specified application depth of Hp : Hence, the pesticide application rate [M/L3/T] can be expressed as MðtÞ ¼ Ms;t =Hp

MSTG;t ¼ ð1 2 Dtkc;t ÞMSTG;t21 þ Mc;t 2 Mcs;t

in which ht is the pesticide wash-off ratio for the plant canopy zone and it is given by Pcs;t þ Ics;t St21 þ Pc;t þ Ic;t 2 Ec;t

ðA9Þ

ðA4Þ

Assuming that the mass of pesticide washed off the plant canopy zone is proportioned to the amount of the throughfall water from the plant canopy to the soil surface, we have   Mcs;t ¼ ht ð1 2 Dtkc;t ÞMSTG;t21 þ Mc;t ðA5Þ

ht ¼

37

ðA6Þ

where St21 is the water storage in the plant canopy zone at time t 2 1; Pc;t is the portion of rainfall on the plant canopy at time t; Ic;t is the portion of overcanopy irrigation on the plant canopy at time t; and Ec;t is the evaporation from the plant canopy at time t: The aforementioned method is equivalent to calculating the pesticide concentration in the plant canopy zone at time t based on the total mass of pesticide ½ð1 2 Dtkc;t ÞMSTG;t21 þ Mc;t  and the total amount of water ðSt21 þ Pc;t þ Ic;t 2 Ec;t Þ: Then the quantity of pesticide washed off the plant canopy is equal to the product of the resulting concentration and the corresponding throughfall water ðPcs;t þ Ics;t Þ: Thus, the storage of pesticide in the plant canopy zone at time t can be written as

Appendix B. Pesticide runoff, erosion, volatilization, and root uptake Following Mullins et al. (1993), we use SCS-CN method, MUSLE, and boundary layer approach to estimate pesticide runoff, erosion, and volatilization, respectively, in our semidiscrete pesticide transport model. B.1. Pesticide runoff The SCS-CN method can be applied to estimate rainfall-induced surface runoff (Singh, 1992) Q¼

ðP 2 0:2SÞ2 P þ 0:8S

ðB1Þ

in which S ¼ ð2540=CNÞ 2 25:4

ðB2Þ

where S is the potential retention (cm); CN, the runoff curve number; P; the precipitation (cm); Q is the runoff (cm). The runoff flux of dissolved pesticide Fr [M/L2/T] from the crop field can be calculated by Fr ¼ ar QCl ¼ Pr Cl

ðB3Þ

in which MSTG;t ¼ ð1 2 ht Þbð1 2 Dtkc;t ÞMSTG;t21 þ Mc;t c

ðA7Þ

Hence, the total soil pesticide application on the ground surface at time t is given by Ms;t ¼ Mds;t þ Mcs;t   ¼ Mds;t þ ht ð1 2 Dtkc;t ÞMSTG;t21 þ Mc;t

Pr ¼ a r Q

where ar is the unit conversion factor; Cl ; the dissolved-phase pesticide concentration [M/L3]; Pr is the pesticide runoff rate coefficient [L/T]. Thus, the pesticide runoff rate [1/T] in Eq. (1) is given by rr ¼ Pr =Dz

ðA8Þ

In the case of over-canopy and soil surface application of pesticides, it is assumed that the total pesticide applied on the soil surface is directly distributed in the thin surface zone (the top discretized cell). For soil-incorporated pesticide applications,

ðB4Þ

ðB5Þ

where Dz is the discretized spatial increment [L]. B.2. Pesticide erosion Based on the Universal Soil-Loss Equation (USLE), Williams (1975) proposed the following

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

38

MUSLE to predict eroded sediment yield for individual storms (Singh, 1992; Mullins et al., 1993) Ye ¼ aye ðVQ qp Þ0:56 Ke ðLSÞCc Pc

ðB6Þ

where Ye is the erosion sediment loss (metric ton/day); aye ; the unit conversion factor ðaye ¼ 11:8Þ; VQ ; the daily volume of event runoff (m3); qp ; the peak storm discharge (m3/s); Ke ; the soil erodibility factor; LS, the topographic factor; Cc ; the soil cover factor; Pc is the conservation practice factor. The erosion flux of the adsorbed-phase pesticide Fe [M/L2/T] can be calculated by taking into account the enrichment ratio for the organic matter as follows Fe ¼

ae Ye rom Cs aYr K C ¼ e e om d l ¼ Pe Cl Aw Aw

ðB7Þ

in which aYr K Pe ¼ e e om d Aw

ðB8Þ

where ae is the unit conversion factor; rom ; the enrichment ratio for organic matter [M/M]; Cs ; the adsorbed-phase pesticide concentration [M/M]; Kd ; the distribution coefficient [L3/M]; Aw ; the watershed area [L2]; Pe is the pesticide erosion rate coefficient [L/T].Thus, the pesticide erosion rate [1/T] in Eq. (1) is given by re ¼ Pe =Dz

ðB9Þ

the vapor-phase pesticide concentration at the soil surface [M/L3], Cga is the vapor-phase pesticide concentration at the top of the air boundary layer [M/L3], and d is the thickness of the stagnant air boundary layer [L]. By assuming Cga ¼ 0 and substituting Cg ð0; tÞ ¼ KH Cl ð0; tÞ; the volatilization flux can be rewritten as Fv ¼ 2

Da KH Cl ð0; tÞ ¼ 2Pv Cl ð0; tÞ d

ðB12Þ

in which Pv ¼ Da KH =d

ðB13Þ

B.4. Plant root uptake By assuming passive root uptake, the following approach, used by Hantush and Marin˜o (1996), is applied for simulating plant uptake of pesticides in the developed model. Boesten and van der Linden (1991) expressed the rate of pesticide uptake by plants in the root zone [1/T] as ru ¼ FS

ðB14Þ

in which S¼a

ETp ½1 2 expð20:6IÞ Lr

ðB15Þ

where arom is the unit conversion factor.

where F is the transpiration stream concentration factor, S is the actual rate of water uptake by the crop [1/T], a is the reduction factor, ETp is the potential rate of evapotranspiration [L/T], Lr is the depth of the crop root zone [L], and I is the leaf area index [L2/L2]. The transpiration stream concentration factor F can be approximated by (Briggs et al., 1982)

B.3. Volatilization

F ¼ 0:784 exp½2ðlog Kow 2 1:78Þ2 =2:44

The enrichment ratio can be estimated by using an empirical approach (Mockus, 1972; Mullins et al., 1993) lnðrom Þ ¼ 2 2 0:2 lnðarom Ye =Aw Þ

ðB10Þ

By assuming that a vapor-phase pesticide diffuses through a stagnant air boundary layer from the soil surface to the overlying atmosphere, Jury et al. (1983, 1991) expressed the volatilization flux Fv [M/L2/T] as Cg ð0; tÞ 2 Cga Fv ¼ 2Da d

ðB11Þ

where Da is the vapor-phase pesticide diffusion coefficient in the air layer [L2/T], Cg ð0; tÞ is

ðB16Þ

where Kow is the octanol –water distribution coefficient [L3/M].

References Barbash, J.E., 1995. Pesticide in ground water: current understanding of distribution and major influences, US Geological Survey Fact Sheet FS-244-95, 4 pp, Available at: http://water. wr.usgs.gov/pnsp/gw/

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40 Bear, J., 1972. Dynamics of Fluids in Porous Media. New York. Beltman, W.H.J., Boesten, J.J.T.I., van der Zee, S.E.A.T.M., 1995. Analytical modeling of pesticide transport from the soil surface to a drinking water well. J. Hydrol. 169, 209– 228. Beltman, W.H.J., Boesten, J.J.T.I., van der Zee, S.E.A.T.M., Quist, J.J., 1996. Analytical modeling of effects of application frequency on pesticide concentrations in wells. Ground Water 34, 470 –479. Bicknell, B.R., Imhoff, J.C., Kittle, J.L. Jr., Donigian, A.S. Jr., Johanson, R.C., 1996. Hydrological Simulation Program— FORTRAN (HSPF), User’s Manual for Release 11, US Environmental Protection Agency, Athens, GA. Boesten, J.J.T.I., van der Linden, A.M.A., 1991. Modeling the influence of sorption and transformation on pesticide leaching and persistence. J. Environ. Qlty 20, 425–435. Briggs, G.G., Bromilow, R.H., Evans, A.A., 1982. Relationships between lipophilicity and root uptake and translocation of nonionised chemicals by barley. Pesticide Sci. 13, 495–504. Carsel, R.F., Smith, C.N., Mulkey, L.A., Dean, J.D., Jowise, P., 1984. User’s Manual for the Pesticide Root Zone Model (PRZM): Release 1, US Environmental Protection Agency, Athens, GA, EPA-600/3-84-109. Carsel, R.F., Nixon, W.B., Ballantine, L.G., 1986. Comparison of pesticide root zone model predictions with observed concentrations for the tobacco pesticide metalaxyl in unsaturated zone soils. Environ. Toxicol. Chem. 5, 345–353. Carsel, R.F., Imhoff, J.C., Hummel, P.R., Cheplick, J.M., Donigian, A.S. Jr., 1998. PRZM3: A Model for Predicting Pesticide and Nitrogen Fate in the Crop Root and Unsaturated Soil Zones: Users Manual for Release 3.0, National Exposure Research Laboratory, US Environmental Protection Agency, Athens, GA. Chu, X., Basagaoglu, H., Marin˜o, M.A., Volker, R.E., 2000. Aldicarb transport in subsurface environment: Comparison of models. J. Environ. Engng ASCE 126 (2), 121–129. Dean, J.D., Huyakorn, P.S., Donigian, A.S., Voss, K.A., Schanz, R.W., Meeks, Y.T., Carsel, R.F., 1989. Risk of Unsaturated/ Saturated Transport and Transformation of Chemical Concentrations (RUSTIC), Theory and Code Verification, vol. 1. US Environmental Protection Agency, Athens, GA, EPA-600/3-89/ 048a. Domagalski, J.L., 1997. Results of a prototype surface water network design for pesticides developed for the San Joaquin River Basin, California. J. Hydrol. 192, 33–50. Domagalski, J.L., Dubrovsky, N.M., Kratzer, C.R., 1997. Pesticides in the San Joaquin River, California: inputs from dormant sprayed orchards. J. Environ. Qlty 26, 454–465. Dubrovsky, N.M., Kratzer, C.R., Brown, L.R., Gronberg, J.M., Burow, K.R., 1998. Water quality in the San Joaquin-Tulare Basins, California, 1992 – 95, US Geological Survey Circular 1159, Available at http://water.usgs.gov/pubs/ circ1159. Ferrari, C.A., McElhiney, M.A., 1997. Soil survey of Stanislaus County, California, Western Part, US Department of Agriculture, Natural Resources Conservation Service, Available at http://www.ca.nrcs.usda.gov/mlra/wstan/main.html

39

Friedly, J.C., 1972. Dynamic Behavior of Processes, Prentice-Hall, Englewood Cliffs, NJ. Guymon, G.L., 1970. A finite element solution of the onedimensional diffusion –convection equation. Water Resour. Res. 6 (1), 204–210. Guymon, G.L., Scott, V.H., Herrmann, L.R., 1970. A general numerical solution of the two-dimensional diffusion–convection equation by the finite element method. Water Resour. Res. 6 (6), 1611–1617. Hantush, M.M., Marin˜o, M.A., 1996. An analytical model for the assessment for pesticide exposure levels in soils and groundwater. Environ. Modeling Assess. 1 (4), 263 –276. Hantush, M.M., Marin˜o, M.A., Islam, M.R., 2000. Models for leaching of pesticides in soils and groundwater. J. Hydrol. 227, 66 –83. Hantush, M.M., Govindaraju, R.S., Marin˜o, M.A., Zhang, Z., 2002. Screening model for volatile pollutants in dual porosity soils. J. Hydrol. 260, 58–74. Hwang, J.C., Cho, W.C., Yeh, G.T., 1984. An eigenvalue solution continuous in time to the spatially discretized solute transport equation in steady groundwater flow. Water Resour. Res. 20 (11), 1725–1732. Jones, R.L., Black, G.W., Estes, T.L., 1986. Comparison of computer model predictions with unsaturated zone field data for aldicarb and aldoxycarb. Environ. Toxicol. Chem. 5, 1027–1037. Jury, W.A., 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, 2465–2474. Jury, W.A., Spencer, W.F., Farmer, W.J., 1983. Behavior assessment model for trace organics in soil. I. Model description. J. Environ. Qlty 12 (4), 558–564. Jury, W.A., Focht, D.D., Farmer, W.J., 1987. Evaluation of pesticide groundwater pollution potential from standard indices of soil-chemical adsorption and biodegradation. J. Environ. Qlty 16 (4), 422–428. Jury, W.A., Gardner, W.R., Gardner, W.H., 1991. Soil Physics, fifth ed., Wiley, New York. Kamra, S.K., Singh, S.R., Rao, K.V.G.K., van Genuchten, M.Th., 1991a. A semidiscrete model for water and solute movement in tile-drained soils. 1. Governing equations and solution. Water Resour. Res. 27 (9), 2439–2447. Kamra, S.K., Singh, S.R., Rao, K.V.G.K., van Genuchten, M.Th., 1991b. A semidiscrete model for water and solute movement in tile-drained soils. 2. Field validation and applications. Water Resour. Res. 27 (9), 2449–2456. Kidd, H., James, D.R., 1991. The Agrochemicals Handbook, third ed., Royal Society of Chemistry Information Services, Cambridge. Knisel, W.G. (Ed.), 1980. CREAMS: A field-Scale Model for Chemicals, Runoff, and Erosion from Agricultural Management Systems, US Department of Agriculture, Science and Education Administration, Conservation Research Report No. 26, USDA, Washington, 643 pp. Knisel, W.G., Leonard, R.A., Davis, F.M., Nicks, A.D., 1993. GLEAMS Version 2.10. Part III. User Manual, Available at

40

X. Chu, M.A. Marin˜o / Journal of Hydrology 285 (2004) 19–40

http://www.wcc.nrcs.usda.gov/water/quality/common/h2oqual. html Kuivila, K.M., Foe, C.G., 1995. Concentrations, transport and biological effects of dormant spray pesticides in the San Francisco Estuary, California. Environ. Toxicol. Chem. 14 (7), 1141–1150. Larson, S.J., Gilliom, R.J., Capel, P.D., 1999. Pesticides in streams of the United States: initial results from the National Water Quality Assessment program, USGS, Water-Resources Investigations Report 98-4222, Sacramento, California,. Leonard, R.A., Knisel, W.G., Still, D.A., 1987. GLEAMS: groundwater loading effects of agricultural management systems. Trans. ASAE 30, 1403–1418. Le Veque, R.J., 1992. Numerical Methods for Conservation Laws, second ed., Birkhauser, Basel. Melancon, S.M., Pollard, J.E., Herm, S., 1986. Evaluation of SESOIL, PRZM, and PESTAN in a laboratory column leaching experiment. Environ. Toxicol. Chem. 5, 865 –878. Mockus, V., 1972. Estimation of direct surface runoff from storm rainfall. In National Engineering Handbook. Section IV, Hydrology. US Soil Conservation Report. NEH-Notice 4-102. Mullins, J.A., Carsel, R.F., Scarbrough, J.E., Ivery, A.M., 1993. PRZM-2, A Model for Predicting Pesticide Fate in the Crop Root and Unsaturated Soil Zones: Users Manual for Release 2.0, US Environmental Protection Agency. Nalluswami, M., Longenbaugh, R.A., Sunada, D.K., 1972. Finite element method for the hydrodynamic dispersion equation with mixed partial derivatives. Water Resour. Res. 8 (5), 1247–1250. Panshin, S.Y., Dubrovsky, N.M., Gronberg, J.M., Domagalski, J.L., 1998. Occurrence and distribution of dissolved pesticides in the San Joaquin River Basin, California, US Geological Survey, Water-Resources Investigations Report 98-4032, National Water-Quality Assessment, Sacramento. Parrish, R.S., Smith, C.N., Fong, F.K., 1992. Test of the pesticide root zone model and the aggregate model for transport and transformation of aldicarb, metolachlor, and bromide. J. Environ. Qlty 21, 685 –697. Pease, W.S., Albright, D., DeRoos, C., Gottsman, L., Kyle, A.D., Morello-Frosch, R., Robinson, J.C., 1995. Pesticide Contamination of Groundwater in California, University of California, Berkeley. Pennell, K.D., Hornsby, A.G., Jessup, R.E., Rao, P.S.C., 1990. Evaluation of five simulation models for predicting aldicarb and bromide behavior under field conditions. Water Resour. Res. 26, 2679–2693.

Simunek, J., Sejna, M., Van Genuchten, M.Th., 1998. HYDRUS1D, Simulating the One-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Version 2.0, US Salinity Laboratory, USDA. Singh, V.P., 1992. Elementary Hydrology, Prentice Hall, Englewood Cliffs, NJ. Skaggs, T.H., Barry, D.A., 1996. Sensitivity methods for timecontinuous, spatially discrete groundwater contaminant transport models. Water Resour. Res. 32 (8), 2409–2420. Spitz, K., Moreno, J., 1996. A Practical Guide to Groundwater and Solute Transport Modeling, Wiley, New York. Spurlock, F., 1998. Evaluation of current simulation models to predict pesticide movement to ground and surface water under California conditions, Study #177, Department of Pesticide Regulation, Sacramento, CA. Umari, A.M.J., Gorelick, S.M., 1986. The problem of complex eigensystems in the semianalytical solution for advancement of time in solute transport simulations: a new method using real arithmetic. Water Resour. Res. 22 (7), 1149–1154. U.S.G.S., 1998. Pesticide in surface and ground water of the United States: summary of results of the National Water Quality Assessment Program (NAWQA), Pesticides National Synthesis Project, http://ca.water.usgs.gov/pnsp/allsum/index.html Van Genuchten, M.Th., Davidson, J.M., Wierenga, P.J., 1974. An evaluation of kinetic and equilibrium equations for the prediction of pesticide movement through porous media. Soil Sci. Soc. Am. Proc. 38, 29–34. Varshney, P., Tim, U.S., Anderson, C.E., 1993. Risk-based evaluation of ground-water contamination by agricultural pesticides. Ground Water 31, 356–362. Wagenet, R.J., Hutson, J.L., 1989. LEACHM: Leaching Estimation and Chemistry Model, Continuum Water Resources Institute, vol. 2. Version 2, Cornell University, Ithaca, NY. Williams, J.R., 1975. Sediment yield prediction with universal equation using runoff energy factor, In Present and Prospective Technology for Predicting Sediment Yields and Sources. ARS-S-40, US Department of Agriculture, Washington, DC, pp. 244–252. Willis, R., 1979. A planning model for the management of groundwater quality. Water Resour. Res. 15 (6), 1305–1312. Young, R.A., Onstad, C.A., Bosch, D.D., Anderson, W.P., 1994. Agricultural Non-Point Source Pollution Model (AGNPS), Version 4.03, User’s Guide, USDA-ARS, Available at http:// www.wcc.nrcs.usda.gov/water/quality/common/h2oqual.html