Environmental Modelling & Software 22 (2007) 1316e1327 www.elsevier.com/locate/envsoft
IPTM-CS: A windows-based integrated pesticide transport model for a canopyesoil system Xuefeng Chu a,*, Miguel A. Marin˜o b,c a Annis Water Resources Institute, Grand Valley State University, MI 49441, USA Department of Land, Air, and Water Resources, University of California, Davis, CA 95616, USA c Department of Civil and Environmental Engineering, University of California, Davis, CA 95616, USA b
Received 21 May 2005; received in revised form 13 May 2006; accepted 25 August 2006 Available online 27 October 2006
Abstract A Windows-based integrated pesticide transport modeling system (IPTM-CS) has been developed for simulating three-phase (dissolved, adsorbed, and vapor phases) pesticide environmental fate in a canopyesoil system (or in the vadose zone alone). The modeling system integrates pre-processing of data (parameters estimation and input), model run, and post-processing of simulation results (summary tables, Excel spreadsheets, and graphs) in a user-friendly Windows interface. To facilitate parameter estimation, an extensive data supporting system that includes convenient parameter calculators and databases has also been developed and incorporated in the IPTM-CS. The data supporting system also provides links to a number of web-based databases maintained by government agencies and institutions. A time-continuous and space-discrete (TC-SD) method is employed to solve the transport problem and five different solution schemes of varying accuracy and features have been incorporated in the interfaced IPTM-CS. Finally, testing of the IPTM-CS at a field site in Chico, California is presented. Comparison of the simulated and observed diazinon concentrations is conducted for three selected depths (0.15, 0.25, and 0.5 m) and the results are further evaluated by using three methods (linear regression, normalized objective function, and modeling efficiency). The quantitative evaluation indicates that the IPTM-CS yields fairly good simulations. Ó 2006 Elsevier Ltd. All rights reserved. Keywords: IPTM-CS; Transport modeling; Pesticide contamination; Vadose zone; Canopyesoil system; Semidiscrete method
Software availability
1. Introduction
Name of software: Integrated Pesticide Transport Model for a CanopyeSoil System (IPTM-CS) Developer: Xuefeng Chu Contact information: Annis Water Resources Institute, Grand Valley State University, 740 West Shoreline Drive, Muskegon, MI 49441, USA. Tel.: þ1 616 331 3987; fax: þ1 616 331 3864. E-mail:
[email protected] Year first available: 2004 Operating system: Windows 98, NT, 2000, or XP Program language: VB, FORTRAN Program size: 5.63 MB Availability and cost: Free (available via E-mail)
Various models have been developed for simulating pesticide transport and transformation in the subsurface environment. According to solution techniques, they can be categorized as two major types: analytical models (e.g., Jury et al., 1983, 1987; Jury and Gruber, 1989; Beltman et al., 1995; Hantush and Marin˜o, 1996; Chu et al., 2000) and numerical models (e.g., Knisel, 1980; Carsel et al., 1984; Leonard et al., 1987; Dean et al., 1989; Wagenet and Hutson, 1989; Knisel, 1993; Mullins et al., 1993; Knisel and Davis, 2000; Ahuja et al., 2000; Carsel et al., 2003). Analytical models generally require few input data and provide continuous and exact closed-form solutions. They are efficient in computations and easy to use. However, in order to derive analytical solutions, more assumptions are commonly introduced for simplification of the initial and boundary conditions, flow
* Corresponding author. Tel.: þ1 616 331 3987; fax: þ1 616 331 3864. E-mail address:
[email protected] (X. Chu). 1364-8152/$ - see front matter Ó 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2006.08.006
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
conditions, porous media, as well as physical and (bio)chemical processes of the simulated pesticides. Thus, analytical models are often used for screening purposes only. In contrast, numerical models are capable of dealing with complex flow and contaminant transport problems (such as arbitrary initial distribution, various boundary conditions, heterogeneous media, unsteady flow, and time-variant pesticide property parameters). The representative numerical models that are widely used for simulating pesticide transport in the vadose zone include: PRZM (Pesticide Root Zone Model, Carsel et al., 1984, 2003; Mullins et al., 1993), 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, 1993; Knisel and Davis, 2000), LEACHM (Leaching Estimation and Chemistry Model, Wagenet and Hutson, 1989), and RZWQM (Root Zone Water Quality Model, Ahuja et al., 2000). PRZM has also been incorporated in the FOCUS models of the EU (FOCUS: Forum for the Co-ordination of Pesticide Fate Models and their Use). Additionally, three other FOCUS models: MACRO (A Model of Water Flow and Solute Transport in Macroporous Soil, Jarvis et al., 1997; Larsbo and Jarvis, 2003), PEARL (Pesticide Emission Assessment at Regional and Local Scales, Tiktak et al., 2000), and PELMO (Pesticide Leaching Model, Klein, 1995) have also been widely used for assessing pesticide leaching potentials in many European countries. Most recently, a great effort has been made in the GIMMI project, funded by the EU IST program, towards integrating involvement of data providers, scientists, service providers, and final users and providing an inter-operable network of Geographic Information (GI) and GI-based web-services in pesticide modeling and environmental impact assessment (Denzer et al., 2005). Similarly, Henriksen et al. (in press) examined application of Bayesian networks (Bns) as an environmental decision support tool for public participation modeling in the management of groundwater contamination induced by pesticides. Use of pesticide transport models generally requires certain modeling experience and computer programming skills, especially if the models do not provide users with any Windowsbased interface. Additionally, estimation of model parameters is often a tedious and time-consuming work. Hence, a stateof-the-art, Windows-based pesticide transport model, equipped with a database for parameter estimation can be very useful. The objective of this study is to develop such an integrated pesticide transport model for a canopyesoil system (IPTMCS) that possesses a flexible model structure, a user-friendly Windows interface, and an extensive data supporting system. The first feature is achieved by using a semidiscrete solution method that integrates analytical and numerical techniques (Chu and Marin˜o, 2004). The model thus can be used as either a screening-level analytical model, such as a three-zone (surface zone, crop root zone, and deep vadose zone) analytical model (e.g., a two-zone analytical model for simulating pesticide transport in the crop root zone and intermediate vadose
1317
zone developed by Hantush and Marin˜o, 1996), or a numerical model that accounts for detailed distribution and variability in pesticide concentrations along the soil profile (e.g., many PRZM applications). To enhance the applicability of the model to a variety of pesticide transport problems and improve the accuracy and stability of the simulations, five different solution schemes are incorporated in the modeling system. The second and third features of the modeling system (Windowsbased interface and data supporting system) especially facilitate file management, data preparation, model execution, post-processing, and parameter estimation. Unlike other comprehensive modeling systems that cover various water quality problems, the IPTM-CS focuses on pesticide transport modeling and information only. Development of this Windows-based IPTM-CS is primarily based on our previous work supported by the USEPA Center for Ecological Health Research (CEHR) at the University of California, Davis. Chu and Marin˜o (2004) described the theoretical details on the model development and the time-continuous and space-discrete solution method. This paper focuses on development of the Windows-based modeling software. For continuity, however, some theoretical background regarding the model development is briefly summarized in the following section. 2. Model development e mathematical expressions and solution methods 2.1. Pesticide fate and transport in the canopy and soil system The plant canopy zone plays an important role in the fate and transport of pesticides by altering their environmental pathways and spatial and temporal distributions for either over-canopy or under-canopy pesticide applications. For the latter case, the pesticide fate is indirectly affected by the plant canopy by changing hydrologic variables (water flow and distributions). If ignoring the portion of pesticide that is directly volatized into the atmosphere during applications, the total application of pesticide sprayed over the canopy (M ) can be partitioned into canopy application (Mc, the portion of pesticide applied on the plant canopy) and direct soil application (Mds, the portion of pesticide directly applied to the soil surface). The canopy-applied pesticide (Mc) can be further washed off the plant canopy zone to the soil surface during rainfall or over-canopy irrigation (Mcs, secondary pesticide application to the soil surface). The pesticide in the canopy zone is also subject to decay and volatilization. The processes that pesticide residues undergo in the canopy zone are shown in Fig. 1 and the overall mass balance of pesticide in the plant canopy zone can be expressed as, dMSTG ðtÞ ¼ mc ðtÞ mcs ðtÞ kc ðtÞMSTG ðtÞ dt
ð1Þ
where MSTG(t) is the pesticide mass storage in the plant canopy zone [M/L2], mc(t) is the canopy pesticide application rate [M/L2/T], mcs(t) is the canopyesoil pesticide application
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
1318
Atmosphere
rq ¼ ð1 KH Þ Total Over-Canopy Pesticide Application M Decay & Volatilization Mc Plant Canopy Zone
MSTG Mcs
Volatilization Infiltration
Mds
Runoff & Erosion
Soils
Fig. 1. Pesticide fate and transport in the plant canopy zone.
rate [M/L2/T], and kc(t) is the first-order pesticide decay/volatilization rate in the plant canopy zone [1/T]. The canopy zone model is solved and the resulting pesticide transfer from the canopy to the soil surface is then incorporated in the pesticide loading term of the following soil model. By taking into account pesticide fate in the canopy zone, the model is able to deal with various pesticide application methods, including over-canopy spray, under-canopy soil surface spray, and any combined over- and under-canopy applications. For soil-incorporated applications of pesticide, the applied pesticide can be directly added to the involved spatial cells within the specified application depth. By assuming linear equilibrium sorption, linear equilibrium liquidevapor partitioning, and first-order decay, the secondorder partial differential equation (PDE) governing one-dimensional three-phase pesticide transport in the vadose zone can be expressed as (Chu and Marin˜o, 2004), vC v v v g vðKH CÞ l vC qR ¼ E þ E ðqCÞ rC þ M ð2Þ vt vz vz vz vz vz in which, R ¼ 1 þ ðrKd þ aKH Þ=q
ð3Þ
Eg ¼ a10=3 =n2 Dag
ð4Þ
El ¼ q10=3 =n2 Dwl þ aL jqj
ð5Þ
r ¼ rr þ re þ ru þ rd þ rq
ð6Þ
r r ¼ ar re ¼
Q Ls
ae Ye rom Kd Ls A
ð7Þ
ð8Þ
ru ¼ FSc
ð9Þ
rd ¼ qRks
ð10Þ
vq vt
ð11Þ
where C is the concentration of the dissolved-phase pesticide [M/L3]; q is the water flux [L/T]; M is the pesticide loading term [M/L3/T]; R is the retardation factor; r is the bulk density [M/L3]; Kd is the distribution coefficient [L3/M]; KH is the dimensionless Henry’s law constant; q is the volumetric water content [L3/L3]; a is the volumetric air content [L3/L3]; n is the porosity [L3/L3]; Dag is the binary diffusion coefficient of the vapor-phase pesticide in free air [L2/T]; Dw l is the binary diffusion coefficient of the dissolved-phase pesticide in water [L2/T]; aL is the longitudinal dispersivity [L]; ks is the first-order degradation rate [1/T]; rr is the pesticide runoff rate [1/T]; re is the pesticide erosion rate [1/T]; ru is the pesticide root uptake rate [1/T]; Ls is the thickness of the surface zone [L]; Q is the runoff [L/T]; Ye is the eroded sediment yield [M/T]; rom is the enrichment ratio for organic matter [M/M]; A is the watershed area [L2]; F is the transpiration stream concentration factor (Briggs et al., 1982; Boesten and van der Linden, 1991); Sc is the actual rate of water uptake by the crop [1/T]; and ar and ae are unit conversion factors. Note that in the IPTM-CS, the SCS-CN method (USDA, 1986) is used to estimate surface runoff and three modified universal soil loss equations (MUSLE, MUSS, and MUST) (Williams, 1975, 1995) are utilized to estimate the sediment yield. The governing Eq. (2) is subject to the following initial and boundary conditions: Cðz; tÞ ¼ C0 ðzÞ ðt ¼ 0Þ vC Da þ q þ KH E KH þ E vz d Da C ¼ qCin þ Cag ðt > 0; z ¼ 0Þ d
g
ð12Þ
l
vC þ qC ¼ qCout ðt > 0; z ¼ Ls Þ Eg KH þ El vz
ð13Þ
ð14Þ
where Da is the vapor-phase pesticide diffusion coefficient in the air layer [L2/T]; Cag is the vapor-phase pesticide concentration at the top of the air boundary layer [M/L3]; Cin is the dissolved-phase pesticide concentration of the inflow at the upper boundary (z ¼ 0) [M/L3]; Cout is the dissolved-phase pesticide concentration of the outflow at the lower boundary (z ¼ Ls) [M/L3]; and d is the thickness of the stagnant air boundary layer [L]. For the upper boundary condition expressed in Eq. (13), it is assumed that a vapor-phase pesticide diffuses through a stagnant air boundary layer from the soil surface to the overlying atmosphere (Jury et al., 1991). 2.2. Semidiscrete solution method A time-continuous and space-discrete semidiscrete solution method is used in the IPTM-CS to solve the secondorder PDE [Eq. (2)] together with the initial and boundary conditions [Eqs. (12)e(14)]. Thus, only the space domain
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
1319
Fig. 2. Main interface of IPTM-CS.
that includes the surface zone, plant root zone and deep vadose zone is discretized. To enhance the model applicability to a variety of problems, different schemes of varying numerical properties (e.g., accuracy and stability) are incorporated in the IPTM-CS modeling system, which can be selected to approximate the spatial derivatives as follows for the advection term: (1) Two-point, first-order accurate, upwind finite difference: vðqCÞ qi Ci qi1 Ci1 z vz i Dz
ð15Þ
(2) Two-point, central differencing of second-order accuracy: vðqCÞ 1 z ðqCÞiþ1 ðqCÞi1 vz i 2Dz
ð16Þ
(3) Four-point, upwind-biased scheme e QUICK (Quadratic Upstream Interpolation for Convective Kinematics) (Leonard, 1979): vðqCÞ 1 ðqCÞi2 7ðqCÞi1 þ3ðqCÞi þ3ðqCÞiþ1 z vz 8Dz i
Fig. 3. Interface of the ParametereEstimation menu.
ð17Þ
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
1320
Fig. 4. Example of the Windows interface for pesticide properties parameters in the data supporting system.
(4) Four-point, third-order upwind-biased scheme: vðqCÞ 1 z ðqCÞi2 6ðqCÞi1 þ3ðqCÞi þ2ðqCÞiþ1 vz i 6Dz
Application of the semidiscrete solution method to the original PDE [Eq. (2)] results in a system of ordinary differential equations (ODEs) in the following matrix form: ð18Þ
(5) van Leer flux limiter (van Leer, 1974): vðqCÞ ðqCÞi ðqCÞi1 zwvl vz i Dz
ð19Þ
C_ ¼ AC þ M
ð23Þ
where A is the coefficient matrix, C is the concentration vector, C_ is the derivative vector of the concentration, and M is the pesticide loading vector. The solution for instantaneous and continuous pesticide loadings can be, respectively, given by (Chu and Marin˜o, 2004):
in which, 4ðri Þ 4ðri1 Þ wvl ¼ 1 þ 2 2ri1 ri þ jri j 1 þ jri j
ð21Þ
ðqCÞiþ1 ðqCÞi ðqCÞi ðqCÞi1
ð22Þ
4ðri Þ ¼
ri ¼
ð20Þ
Chu and Marin˜o (2006) examined the accuracy, applicability, and limitations of these five numerical schemes in detail and compared their performances with each other and also with standard numerical and analytical counterparts in a series of testing problems.
CðtÞ ¼ eAðtt0;i Þ ½Cðt0;i Þ þ m
ð24Þ
CðtÞ ¼ eAðtt0;i Þ Cðt0;i Þ þ A1 eAðtt0;i Þ I m
ð25Þ
The reference provided detailed derivations and a review of applications of the semidiscrete approach to contaminant transport modeling in porous media. Evaluating the matrix exponential is an essential part of the model, which further depends on the properties of the coefficient matrix A. Generally, any scheme that tends to weaken or eliminate the diagonal dominance of the matrix A may lead to instability in the solution. This issue and computation methods for the matrix exponential have been discussed in depth by Moler and Van Loan (1978) and Ogata (1987).
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
1321
Fig. 5. Example of pesticide parameter calculators in the data supporting system (distribution coefficient calculator).
3. Software development e windows-based IPTM-CS To facilitate pesticide transport modeling for the canopye soil system described in the preceding section, a Windowsbased integrated pesticide transport modeling system, IPTM-CS, has been developed. The IPTM-CS main interface (Fig. 2) is designed according to three fundamental modeling procedures: data preparation, model execution, and post-processing, which are covered by menus Data, Model, and Output, respectively. The menu Data covers eight types of input data/ parameters associated with the period of simulation, space domain and discretization, soil hydraulic properties, pesticide properties and applications, plant canopy features, climatic conditions, hydrologic characteristics and irrigation, and solution methods. The IPTM-CS supports flexible data input methods. All time- and/or space-varying data, such as rainfall data, can be either inputted manually or imported directly from a text-formatted file or a Microsoft Excel file. The input data can also be exported and saved to an external file. The five different solution methods, described in the preceding section, have been incorporated in the IPTM-CS for user selection. Generally, the first-order upwind scheme can be the first choice for most problems since it possesses good stability. Yet, it may introduce significant numerical dispersion. The van Leer flux limiter is also a stable numerical scheme and it is especially suitable for advection-dominated transport problems. But it
may take much more computation time. For the central differencing scheme and the two other multi-point, high-order, upwind-biased methods, stability (artificial oscillations) can be an issue for some problems although they yield more accurate solutions. Thus, they have strict application conditions. In the IPTM-CS, irrigation data can be either inputted into the model if actual irrigation data are available or generated by the automatic irrigation simulators in the IPTM-CS. The built-in
Nest #6
Nest #3
Nest #5
Nest #2
3.05m Nest #4
Nest #1 4.57m
Fig. 6. Layout map of Nests no. 1e6.
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
1322
Table 1 Observed diazinon concentrations at the three depths (0.15, 0.25, and 0.50 m) Date
Nest no.
C (mg/L) at 0.15 m
Date
Nest no.
C (mg/L) at 0.25 m
Date
Nest no.
C (mg/L) at 0.50 m
3/4/01 3/6/01 3/8/01 3/11/01 3/23/01 3/6/01 3/8/01 3/11/01 3/14/01 3/20/01
1 1 1 1 1 2 2 2 2 2
0.1828 0.2522 0.1505 0.1486 0.0019 0.2819 0.0751 0.0924 0.0361 0.0060
3/6/01 3/8/01 3/11/01 3/14/01 3/17/01 3/20/01 3/23/01 3/27/01 3/29/01 3/8/01 3/11/01 3/14/01 3/17/01 3/23/01 3/27/01 3/6/01 3/8/01 3/11/01 3/14/01 3/17/01 3/20/01 3/23/01 3/27/01 3/23/01 3/24/01 3/27/01 4/4/01
3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 5 5 5 5 5 5 5 5 6 6 6 6
0.1029 0.0708 0.0510 0.0179 0.0319 0.0052 0.0406 0.0073 0.0114 0.0894 0.0526 0.0861 0.0333 0.0480 0.0458 0.1616 0.0434 0.0492 0.0472 0.0195 0.0152 0.0049 0.0040 0.0224 0.0113 0.0167 0.0120
3/14/01 3/22/01 3/23/01 3/27/01 3/29/01 3/6/01 3/11/01 3/14/01 3/17/01 3/20/01 3/27/01 3/29/01 4/4/01 4/13/01 3/8/01 3/11/01 3/14/01 3/17/01 3/20/01 3/23/01 3/27/01 3/29/01 4/4/01
2 2 2 2 2 3 3 3 3 3 3 3 3 3 5 5 5 5 5 5 5 5 5
0.0546 0.0052 0.0195 0.0048 0.0019 0.0746 0.0276 0.0682 0.0597 0.0044 0.0048 0.0211 0.0135 0.0190 0.0670 0.0554 0.0674 0.0190 0.0199 0.0152 0.0030 0.0123 0.0166
properties, plants, agricultural practices, hydrologic processes, as well as meteorological information. In addition, some calculators have been developed and they are available in the ParametereEstimation menu for computing the corresponding parameters, such as distribution coefficient, Henry’s constant, and first-order decay rate (Fig. 5). The ParametereEstimation menu also provides links to a number of web-based pesticide, soil, plant, and climate databases maintained by government agencies and institutions [e.g., USDA-ARS Pesticide Properties Database, EXTOXNET Pesticide Information Profiles (PIPs), USDA-NRCS National Soil Survey Geographic (SSURGO) Database, State Soil Geographic (STATSGO) Database, and NOAA National Climate Data Center
2.5 2.0
Rainfall (cm)
irrigation simulators facilitate irrigation scheduling according to soil moisture conditions and a variety of management criteria in terms of water use efficiency and water quality control. Both over-canopy and under-canopy irrigation methods and both runoff and no-runoff irrigation management schemes can be implemented to enhance the capability of the model and to satisfy a variety of practical needs for scheduling both timing and quantity of irrigation. The automatic irrigation simulators can be useful for the selection of irrigation methods and identification of appropriate irrigation management strategies. The menu Model includes input data checking and model execution. Through the menu Output, both temporal and spatial distributions of the simulated pesticide concentrations can be shown in tabular and graphical formats. The simulated concentrations can also be exported to Microsoft Excel for further processing. Additionally, this Output menu also provides a number of summary tables, such as water and pesticide tables for the plant canopy zone and the soil surface zone, soil water content table, as well as pesticide runoff and erosion tables. To facilitate parameter estimation, a comprehensive data supporting system has been developed and incorporated in the IPTM-CS. It can be accessed via menu ParametereEstimation (Fig. 3). All parameters/data used in the model are detailed in the data supporting system. The system covers all basic data information (e.g., definition of the model parameters/data and estimation approaches) and includes a number of tables and figures that provide the specific values/ranges of the parameters (Fig. 4) related to pesticide and soil
1.5 1.0 0.5 0.0 2/1/01
3/1/01
4/1/01
5/1/01
Date Fig. 7. Rainfall at Station CHICO.T (Touchtone #08, Chico, CA), 2/1/2001e4/ 30/2001.
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327 Table 2 Soil hydraulic parameters
0.8
Saturated hydraulic conductivity Ks (m/d) Saturated water content qs Residual water content qr Field capacity qFC Wilting point qWP Initial water content q0 Soil water retention parameter n Porosity ns Bulk density r (g/cm3)
0.3
1.061
0.43
0.41
0.078
0.065
0.23 0.117 0.2
0.207 0.095 0.15
1.56
1.89
0.43 1.42
0.41 1.49
Carsel and Parrish (1988) Carsel and Parrish (1988) Carsel and Parrish (1988) Carsel et al. (2003) Carsel et al. (2003)
Carsel and Parrish (1988) Carsel et al. (2003) Carsel et al. (2003)
(NCDC)]. These databases provide much more detailed information on pesticides, soils, plants, and climate. 4. Model testing e application in Chico, California Diazinon, an organophosphate insecticide, is widely used in orchards in the Central Valley of California. Due to intense applications of diazinon especially during rainfall seasons, high concentrations of diazinon, primarily resulting from the storm runoff and erosion, have been frequently detected in the Sacramento and San Joaquin Rivers as well as in their tributaries (Domagalski et al., 2000; Kratzer et al., 2002). Diazinon and chlorpyrifos have been identified as the primary pesticides threatening the quality of surface water in California. Compared with organochlorine pesticides, such as DDT, diazinon has much lower persistence and sorption in soil. Its half-life is about 40 days and the organic carbon partition coefficient (Koc) is 1000 cm3/g (the half-life and Koc of DDT can be as high as 2000 days and 2,000,000 cm3/g, respectively) (Hornsby et al., 1996). Thus, in the subsurface environment, high diazinon concentrations might be limited within shallow soils for most cases (Chu and Marin˜o, 2004). For the aforementioned reasons, we selected diazinon for the model testing. Specifically, a French prune orchard site in Chico, California Table 3 Major parameters on diazinon properties Parameter
Value
References
Distribution coefficient Kd (cm3/g) Henry’s constant KH (dimensionless) First-order decay rate ks (1/day) Diffusion coefficient in free air Dag (m2/day) Diffusion coefficient in water Dwl (m2/day) Air boundary layer thickness d (m) Log of octanolewater partitioning coefficient Log Kow (cm3/g)
0.65e3.23 5.0 105 0.08e0.10 0.43
Carsel et al. (2003) Carsel et al. (2003) Carsel et al. (2003), c Jury et al. (1991)
0.000043
Jury et al. (1991)
0.005 3.02
Jury et al. (1991) Carsel et al. (2003)
0.7
Concentration (mg/L)
D ¼ 0e0.36 m D ¼ 0.36e1.31 m References
0.05m 0.07m 0.09m 0.11m 0.15m 0.21m 0.25m 0.31m 0.35m
0.6 0.5 0.4 0.3 0.2 0.1 0.0 2/1/01
3/1/01
4/1/01
5/1/01
Date Fig. 8. Temporal distribution of the simulated diazinon concentrations for depths 0.05e0.35 m.
(1012 m2) was selected where diazinon was applied by using a spray rig twice at a rate of 0.52 g/m2 on February 13 and February 28, 2001, respectively. Lysimeters were installed at six locations (Nests 1e6, Fig. 6) and various diazinon samples were collected at different soil depths (0.15, 0.25, and/or 0.5 m) following the second application. The specific sampling dates and the corresponding observed concentrations are listed in Table 1 (the table only includes data used for the model testing). The purpose of this section is to present an application of the Windows-based IPTM-CS software to a real-world problem, in which the simulated diazinon concentrations along the soil profile at various times will be compared with those measured in the field. The soil profile simulated in the modeling is 1.31 m, which consists of loam at the top and sandy loam at the bottom. The thickness of the thin surface zone is 0.02 m and the discretized space increment is 0.02 m for shallow soil (0e0.36 m) and 0.05 m for deeper soil (0.36e1.31 m). The simulation period ranges from February 1, 2001 to April 30, 2001 (89 days). Rainfall data at Station CHICO.T (Touchtone #08, Chico, CA) from February 1 to April 30, 2001 (Fig. 7) are used for the modeling. The plant canopy zone intercepts rainwater 0.12 0.35m 0.44m 0.54m 0.64m 0.74m 0.84m 0.94m 1.04m 1.14m 1.24m
0.10
Concentration (mg/L)
Soil horizons
c, Calibrated.
1323
0.08 0.06 0.04 0.02 0.00 2/1/01
3/1/01
4/1/01
5/1/01
Date Fig. 9. Temporal distribution of the simulated diazinon concentrations for depths 0.35e1.24 m.
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
1324
0.18
Concentration (mg/L) 0.00 0.0
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.16
Depth = 0.25 m
Concentration (mg/L)
0.14 0.2
Depth (m)
0.4 3/3/01 3/4/01 3/5/01 3/7/01 3/10/01 3/15/01 3/20/01 3/25/01 3/30/01
0.6 0.8 1.0
Simulated Observed #3 Observed #4 Observed #5 Observed #6
0.12 0.10 0.08 0.06 0.04 0.02 0.00 2/1/01
3/1/01
1.2
4/1/01
5/1/01
Date
Fig. 10. Distribution of the simulated diazinon concentrations along the soil profile for different times in March 2001.
Fig. 12. Comparison of the simulated and observed diazinon concentrations at depth ¼ 0.25 m.
and the plant interception storage capacity is assumed as 0.1e 0.2 cm. A runoff curve number (CN ) of 72 is chosen for the average antecedent moisture condition (AMC-II). Table 2 shows primary soil hydraulic parameters and Table 3 lists the major parameters concerning diazinon properties. The first-order upwind scheme is selected in this testing problem. It is assumed in the simulation that the initial diazinon concentrations along the entire soil profile are zero. The IPTM-CS provides detailed distributions of diazinon in both time and space. The temporal distributions of the simulated diazinon concentrations for depths 0.05e0.35 m and depth 0.35e1.24 m are shown in Figs. 8 and 9, respectively. Two diazinon peaks corresponding to the two times of application are clearly observed in shallow soils from Fig. 8. After the first application, it takes 5, 7, and 10 days for the concentration peak to reach depths of 0.15, 0.25, and 0.35 m, respectively. For the second application on February 28, however, it takes only 3, 6, and 8 days to reach the three depths because of an immediate rainfall on March 2. It is getting harder to distinguish the two peaks for the soils deeper than 1 m (Fig. 9). After April 1, 2001 (one month following the second application), the diazinon concentrations are smaller than 0.01 mg/L for all soil depths. Fig. 10 shows diazinon variations along the
soil profile for different times. In the first week following the second application (before March 7), diazinon concentrations decrease rapidly in the shallow soil of a depth less than 0.35 m. Higher concentration gradients can be observed in the low permeable shallow soil. The concentration curve for March 30 indicates that diazinon residues are almost evenly distributed along the soil column at a very low level. To test the performance of the model, the temporal distributions of the simulated diazinon concentrations for three depths (0.15, 0.25, and 0.5 m) are compared against the observed ones and the comparison results are shown in Figs. 11e13, respectively. Note that the observation period ranges only from March 2 to April 3. The comparison results for the three data sets (simulated and observed diazinon concentrations for depths of 0.15, 0.25, and 0.5 m) are further evaluated by using three quantitative methods: linear regression, normalized objective function (NOF) (Hession et al., 1994; Kornecki et al., 1999; Fox et al., 2004), and modeling efficiency (EF) (Loague and Green, 1991; Larsbo and Jarvis, 2005). The values of NOF and EF are, respectively, given by: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 1 1X 2 NOF ¼ ð26Þ ðCobs;i Csim;i Þ Cobs n i¼1
0.35 0.08 Simulated Observed #1 Observed #2
0.25 0.20 0.15 0.10 0.05 0.00 2/1/01
Depth = 0.5 m
Depth = 0.15 m
Concentration (mg/L)
Concentration (mg/L)
0.30
Simulated Observed #2 Observed #3 Observed #5
0.06
0.04
0.02
0.00 3/1/01
4/1/01
5/1/01
Date Fig. 11. Comparison of the simulated and observed diazinon concentrations at depth ¼ 0.15 m.
2/1/01
3/1/01
4/1/01
5/1/01
Date Fig. 13. Comparison of the simulated and observed diazinon concentrations at depth ¼ 0.5 m.
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
Depth ¼ 0.15 m Depth ¼ 0.25 m Depth ¼ 0.50 m
Pn
EF ¼ 1 Pi¼1 n
Equation
R2
Csim ¼ 0.5255Cobs þ 0.0414 Csim ¼ 0.8949Cobs þ 0.0136 Csim ¼ 0.6073Cobs þ 0.0110
0.7357 0.6927 0.6327
ðCobs;i Csim;i Þ
i¼1 ðCobs;i
NOF
EF
0.4516 0.5773 0.5234
0.6428 0.5656 0.6317
2
2
Cobs Þ
ð27Þ
where Cobs;i is the observed concentration, Csim;i is the simulated concentration, Cobs is the mean of the observed concentrations, and n is the number of the observed and simulated data set. Note that if all observed concentrations are the same as the simulated ones, the NOF and EF values will be 0 and 1, respectively. The results based on the three evaluation methods are presented in Table 4 and the linear regression between the observed and simulated concentrations for the three data sets is also shown in Fig. 14aec. For the first data set (the observed and simulated concentrations at depth ¼ 0.15 m), the model seems to underestimate the diazinon levels although the highest R2 value (R2 ¼ 0.7357) is achieved (Fig. 14a). Actually, there is a time shift between the observed and simulated curves that increases the differences between the two concentration values, as indicated in Fig. 11. In addition, the lowest NOF (NOF ¼ 0.4516) and highest EF (EF ¼ 0.6428) values are achieved for this data set. Linear regression for the second data set (the observed and simulated concentrations at depth ¼ 0.25 m) results in a slope of 0.8949, intercept of 0.0136, and R2 of 0.6927. This data set results in the highest NOF (NOF ¼ 0.5773) and lowest EF (EF ¼ 0.5656) values. According to the graphical comparisons of the observed and simulated diazinon concentrations (Figs. 11e13) and the quantitative analysis of linear regression, normalized objective function, and modeling efficiency, the IPTM-CS yields fairly good simulations (EF ¼ 0.5656e0.6428, NOF ¼ 0.4516e 0.5773, and R2 ¼ 0.6327e0.7357) and the simulations capture the primary variability and magnitude of the diazinon exposure levels for all three soil depths.
Simulated Concentration (mg/L)
Linear regression
0.30 Depth = 0.15m 0.25 0.20 0.15 0.10 0.05 0.00 0.00
b Simulated Concentration (mg/L)
Data set
a
c
0.05
0.10
0.15
0.20
0.25
0.30
0.18 Depth = 0.25m 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0.00
Simulated Concentration (mg/L)
Table 4 Evaluation of the comparison between the observed and simulated concentrations
1325
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.16
0.18
0.08 Depth = 0.50m
0.06
0.04
0.02
0.00 0.00
0.02
0.04
0.06
0.08
Observed Concentration (mg/L) Fig. 14. Linear regression between the observed and simulated diazinon concentrations for the three selected depths (0.15, 0.25, and 0.5 m).
5. Summary and conclusions A user-friendly, Windows-based integrated pesticide transport modeling software (IPTM-CS) was developed and presented. The IPTM-CS describes the environmental fate of pesticides in the canopyesoil system and simulates a set of transport processes such as advection, diffusion/dispersion, linear equilibrium sorption, linear equilibrium partitioning between vapor and dissolved phases, first-order decay, plant root uptake, volatilization from the soil surface to the overlying atmosphere, as well as pesticide runoff and erosion in the surface zone. The modeling system takes into account various pesticide application methods and incorporates different
semidiscrete solution methods. In particular, the IPTM-CS integrates parameters estimation, pre-processing of data, modeling, and post-processing of simulation results in a userfriendly Windows interface. The incorporated data supporting system provides users with essential information on the model input data/parameters (definitions and estimation methods) and convenient ways (databases, parameter calculators, links to other web-based databases) to estimate all input data/parameters related to soils, pesticides, plants, agricultural practice and management, hydrologic processes, as well as meteorological conditions. Although the IPTM-CS was
1326
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327
developed for a canopyesoil system, it can also be used for the vadose zone alone by properly setting some parameters related to the plant canopy and pesticide applications. For testing purposes, the IPTM-CS was applied to a field site in Chico, California to simulate diazinon fate and transport in the subsurface environment. Three observed and simulated diazinon concentration data sets corresponding to the depths of 0.15, 0.25, 0.5 m were selected and compared. The testing results for the three data sets were further evaluated by using three methods: linear regression, normalized objective function, and modeling efficiency. It can be concluded from the quantitative evaluation that the IPTM-CS captured the primary variability and magnitude of the diazinon exposure levels in the subsurface system. It should be noted that only diazinon was selected for the model testing in this study. It would be definitely of great interest to test the model by simulating the fate and transport of various pesticides of distinct persistence, mobility, sorption, and volatilization properties under dissimilar environmental conditions (e.g., comparison between organophosphate and organochlorine pesticides). Pesticides undergo a series of complex physical, (bio)chemical processes in the plant canopy zone and the vadose zone. In addition to their inherent properties, many other factors (such as macropore flow, Fox et al., 2004; Kamra and Lennartz, 2005; Larsbo and Jarvis, 2005) may also affect their fate and transport in the environment. Inevitably, some processes/mechanisms may not be adequately described in the IPTM-CS, which may limit its applicability to certain cases due to the underlying assumptions (e.g., first-order decay and linear equilibrium sorption were assumed and the lumped SCS-CN model was used for estimating surface runoff). Although these aspects can be improved in the future version of the software, there is always a tradeoff between the model robustness and complexity. Acknowledgements We would like to thank Mr. Nick Lasher for providing valuable field data in Chico, California that have been used for the model test in this study. The development of the Windowsbased IPTM-CS was supported by the Annis Water Resources Institute, Grand Valley State University. The original modeling work was supported by the U.S. Environmental Protection Agency (USEPA) (Grant no. R819658) Center for Ecological Health Research at the University of California, Davis and the UC Toxic Substances Research and Teaching Program. 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. References Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, M.J., Ma, L., 2000. The Root Zone Water Quality Model. Water Resources Publications LLC, Highlands Ranch, CO, 372 pp. 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, 209e228.
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. Qual. 20, 425e435. Briggs, G.G., Bromilow, R.H., Evans, A.A., 1982. Relationships between lipophilicity and root uptake and translocation of non-ionised chemicals by barley. Pestic. Sci. 13, 495e504. Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755e769. 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. U.S. Environmental Protection Agency, Athens, GA (EPA-600/3-84-109). Carsel, R.F., Imhoff, J.C., Hummel, P.R., Cheplick, J.M., Donigian Jr., A.S., 2003. PRZM-3, A Model for Predicting Pesticide and Nitrogen Fate in the Crop Root and Unsaturated Soil Zones: Users Manual for Release 3.12. National Exposure Research Laboratory, USEPA. Chu, X., Marin˜o, M.A., 2004. Semidiscrete pesticide transport modeling and application. J. Hydrol. 285 (1e4), 19e40. Chu, X., Marin˜o, M.A., 2006. Improved compartmental modeling and application to three-phase contaminant transport in unsaturated porous media. J. Environ. Eng. ASCE 132 (2), 211e219. Chu, X., Basagaoglu, H., Marin˜o, M.A., Volker, R.E., 2000. Aldicarb transport in subsurface environment: comparison of models. J. Environ. Eng. ASCE 126 (2), 121e129. 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). In: Theory and Code Verification, vol. 1. U.S. Environmental Protection Agency, Athens, GA (EPA-600/3-89/048a). Denzer, R., Riparbelli, C., Villa, M., Gu¨ttler, R., 2005. GIMMI: geographic information and mathematical models inter-operability. Environ. Model. Softw. 20 (12), 1478e1485. Domagalski, J.L., Knifong, D.L., Dileanis, P.D., Brown, L.R., May, J.T., Connor, V., Alpers, C.N., 2000. Water quality in the Sacramento River basin, California, 1994e98. In: U.S. Geological Survey Circular 1215. USGS, 36 pp., Available from: . Fox, G.A., Malone, R., Sabbagh, G.J., Rojas, K., 2004. Interrelationship of macropores and subsurface drainage for conservative tracer and pesticide transport. J. Environ. Qual. 33, 2281e2289. Hantush, M.M., Marin˜o, M.A., 1996. An analytical model for the assessment for pesticide exposure levels in soils and groundwater. Environ. Model. Assess. 1 (4), 263e276. Henriksen, H.J., Rasmussen, P., Brandt, G., von Bu¨low, D., Jensen, F.V. Public participation modeling using Bayesian networks in management of groundwater contamination. Environ. Model. Softw., in press, doi:10.1016/j.envsoft.2006.01.008. Hession, W.C., Shanholtz, V.O., Mostaghimi, S., Dillaha, T.A., 1994. Uncalibrated performance of the finite element storm hydrograph model. Trans. ASAE 37, 777e783. Hornsby, A.G., Wauchope, R.D., Herner, A.E., 1996. Pesticide Properties in the Environment. Springer, New York. Jarvis, N.J., Hollis, J.M., Nicholls, P.H., Mayr, T., Evans, S.P., 1997. MACRODB: a decision-support tool for assessing pesticide fate and mobility in soils. Environ. Model. Softw. 12 (2e3), 251e265. 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, 2465e2474. Jury, W.A., Spencer, W.F., Farmer, W.J., 1983. Behavior assessment model for trace organics in soil: I. Model description. J. Environ. Qual. 12 (4), 558e564. 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. Qual. 16 (4), 422e428. Jury, W.A., Gardner, W.R., Gardner, W.H., 1991. Soil Physics, fifth ed. John Wiley & Sons, Inc., New York. Kamra, S.K., Lennartz, B., 2005. Quantitative indices to characterize the extent of preferential flow in soils. Environ. Model. Softw. 20, 903e915. Klein, M., 1995. PELMO Pesticide Leaching Model. Version 2.01, User’s Manual. Fraunhofer-Institut, fur Umweltchemie und Okotoxikolgie, D-57392 Schmallenberg.
X. Chu, M.A. Marin˜o / Environmental Modelling & Software 22 (2007) 1316e1327 Knisel, W.G. (Ed.), 1980. CREAMS: A field-scale Model for Chemicals, Runoff, and Erosion from Agricultural Management Systems. Conservation Research Report No. 26. U.S. Department of Agriculture, Science and Education Administration, Washington, 643 pp. Knisel, W.G., 1993. GLEAMS: Groundwater Loading Effects of Agricultural Management Systems. Univ. of Georgia, Coastal Plain Expt. Sta, Bio. and Agri. Engr. Dept., Publ. No. 5, 260 pp. Knisel, W.G., Davis, F.M., 2000. GLEAMS: Groundwater Loading Effects of Agricultural Management Systems. USDA-ARS, Southeast Watershed Research Laboratory, Tifton, Georgia, Version 3.0 Publication No. SEWRLWGK/FMD-050199, revised 081500, 191 pp. Kornecki, T.S., Sabbagh, G.J., Storm, D.E., 1999. Evaluation of runoff, erosion, and phosphorus modeling system: SIMPLE. J. Am. Water Resour. Assoc. 35 (4), 807e820. Kratzer, C.R., Zamora, C., Knifong, D.L., 2002. Diazinon and chlorpyrifos Loads in the San Joaquin River Basin, California, January and February 2000. In: U.S. Geological Survey, Water Resources Investigations Report 02-4103, 38 pp., Available from: . Larsbo, M., Jarvis, N., 2003. MACRO 5.0 A Model of Water Flow and Solute Transport in Macroporous Soil. Technical Description. Swedish University of Agricultural Sciences. Emergo 2003:6 Report, ISSN 1651-7210, ISBN 91-576-6592-3. Larsbo, M., Jarvis, N., 2005. Simulating solute transport in a structured field soil: uncertainty in parameter identification and predictions. J. Environ. Qual. 34, 621e634. Leonard, B.P., 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comp. Meth. Appl. Mech. Eng. 19 (1), 59e98. Leonard, R.A., Knisel, W.G., Still, D.A., 1987. GLEAMS: groundwater loading effects of agricultural management systems. Trans. ASAE 30, 1403e1418.
1327
Loague, K., Green, R.E., 1991. Statistical and graphical methods for evaluating solute transport models: overview and application. J. Contam. Hydrol. 7, 51e73. van Leer, B., 1974. Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second order scheme. J. Comput. Phys. 14, 361e370. Moler, C., Van Loan, C., 1978. Nineteen dubious ways to compute the exponential of a matrix. SIAM Rev. 20 (4), 801e836. 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. U.S. Environmental Protection Agency. Ogata, K., 1987. Discrete-Time Control Systems. Prentice-Hall, Inc., Englewood Cliffs. Tiktak, A., van den Berg, F., Boesten, J.J.T.I., Leistra, M., van der Linden, A.M.A., van Kraalingen, D., 2000. Pesticide Emission Assessment at Regional and Local Scales: User Manual of Pearl Version 1.1. RIVM Report 711401008. In: Alterra Report 28. RIVM, Bilthoven, 142 pp. USDA, 1986. Urban Hydrology for Small Watersheds. Technical Release 55 (TR-55). Natural Resources Conservation Services. Wagenet, R.J., Hutson, J.L., 1989. LEACHM: leaching estimation and chemistry model. Version 2. In: Continuum, vol. 2. Water Resources Inst., Cornell Univ., Ithaca, NY. Williams, J.R., 1975. Sediment yield prediction with universal equation using runoff energy factor. ARS-S-40. In: Present and Prospective Technology for Predicting Sediment Yields and Sources. U.S. Department of Agriculture, Washington, DC, pp. 244e252. Williams, J.R., 1995. The EPIC model. In: Singh, V.P. (Ed.), Computer Models of Watershed Hydrology. Water Resources Publications, Littleton, CO, pp. 909e1000 (Chapter 25).
