Modification of water movement equations in the PRZM3 for simulating pesticides in soil profile

Agricultural Water Management 143 (2014) 38–47

Contents lists available at ScienceDirect

Agricultural Water Management journal homepage: www.elsevier.com/locate/agwat

Modification of water movement equations in the PRZM3 for simulating pesticides in soil profile Masoud Noshadi ∗ , Sajad Jamshidi 1 Department of Water Engineering, College of Agriculture, Shiraz University, Shiraz, Islamic Republic of Iran

a r t i c l e

i n f o

Article history: Received 7 October 2013 Accepted 11 April 2014 Available online 11 July 2014 Keywords: 2,4-D Soil water content Mobile-immobile MC-Cormack

a b s t r a c t In the area of simulating pesticide transportation in soil profile, the majority of models adhere to simplification techniques and due to this, the simulation achieved does not reflect the reality. One such commonly used model in this regard is PRZM3. In the present research, to ensure improved results, in addition to a general software update, modified equations concerning water movement in soil were applied. To achieve this, one of the most exact numerical solutions (MC-Cormack method) was selected for solving water movement equation (Richard’s equation). This equation was then rewritten in the form of mobileimmobile (MIM), and the Shuffled Complex Evaluation (SCE) method for calculating mobile-immobile coefficients was also added to the model. Following model modification, this was used to simulate 2,4-D concentration, and the results were then compared with the results of the main model and measured data (Noshadi et al., 2011) in two different treatments (normal irrigation and deficit irrigation). Considering the statistics, in the normal irrigation treatment for PRZM3, the figure for NRMSE (normalized root mean square error), CRM (coefficient of residual mass) and d (index of agreement) accounted for 0.58, 0.78 and −0.47, respectively while the figures reported in the modified model using MC-Cormack method (PRZM3-MC) were 0.79, 0.28 and −0.04, and in the modified model using MIM form (PRZM3-MC-MIM) they were 0.86, 0.23 and −0.06. Regarding deficit irrigation treatment, for PRZM3, the figure for NRMSE, CRM and d accounted for 0.65, 0.52 and 0.08, respectively while the figures reported in the modified model using PRZM3-MC were 0.77, 0.38 and −0.24 and in PRZM3-MC-MIM they were 0.73, 0.36 and −0.24, respectively. Simulation results reveal that compared to PRZM3, results were more accurate after model modification using PRZM3-MC and PRZM3-MC-MIM. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The phenoxyacetic acid based herbicides constitute a group of compounds used extensively to control the broad-leaf weeds. 2,4-Dichlorophenoxy acetic acid (2,4-D) is commonly preferred because of its low cost and good selectivity. 2,4-D is a colorless solid of melting point = 140.5 ◦ C, vapor pressure = 53 Pa (160 ◦ C), and water solubility = 620 mg L−1 (Worthing and Hance, 1991). Much work has been devoted to investigate the behavior of 2,4-D in soils. Johnson et al. (1995a) found that soil organic content and soil pH are the main determinants of 2,4-D adsorption in soils. Adsorption increases with increasing soil organic content and decreasing soil pH (Johnson et al., 1995a). Degradation rates are determined

∗ Corresponding author. Tel.: +98 7112286130. E-mail addresses: [email protected] (M. Noshadi), [email protected] (S. Jamshidi). 1 Tel.: +98 9376180195. http://dx.doi.org/10.1016/j.agwat.2014.04.011 0378-3774/© 2014 Elsevier B.V. All rights reserved.

by the microbial population, environmental pH, soil moisture, and temperature (Que Hee and Sutherland, 1981; Sandmann et al., 1988; Wilson et al., 1997). In soils, 2,4-D is degraded primarily by microbes. Hemmett and Faust (1969) concluded that the size of the microbial population, the concentration of 2,4-D, and the ratio of the two factors determine 2,4-D degradation rates. Soil conditions that enhance microbial populations (i.e., warm and moist) facilitate 2,4-D degradation (Foster and McKercher, 1973). In addition, 2,4-D has been shown to dissipate more rapidly in soils that were previously treated with 2,4-D, presumably because there was an increase in 2,4-D degrading bacteria after the first application (Oh and Tuovinen, 1991; Smith and Aubin, 1994; Shaw and Burns, 1998). 2,4-D changes form and function with changes in water pH (Que Hee and Sutherland, 1981). In alkaline (high pH; pH > 7) waters, 2,4-D takes an ionized (negatively charged) form that is water-soluble and remains in the water column. Theoretically, in water of a lower pH, 2,4-D remains in a neutral molecular form, increasing its potential for adsorption to organic particles in water, and increasing its persistence (Wang et al., 1994). 2,4-D is most

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

likely to adsorb to suspended particles in muddy waters with a fine silt load (Que Hee and Sutherland, 1981), but little adsorption has been observed in the soil (Halter, 1980). Wilson et al. (1997) found that adequate soil moisture was the most influential parameter affecting degradation rates. Cold, dry soils can hold 2,4-D residues for significantly longer periods (Que Hee and Sutherland, 1981). In relatively dry soils with low bacterial counts, fungi play an increasingly important role in the degradation of 2,4-D (Ou, 1984; Han and New, 1994). Johnson et al. (1995b) found that dissipation rates did not differ significantly between rice field soils and bare ground, suggesting that plants do not play a significant role in eliminating 2,4-D from soils. Moreover, the potential for 2,4-D to volatilize increases with increasing temperature, increasing soil moisture, and decreasing clay and organic matter content in the soil (Helling et al., 1971). Johnson et al. (1995a,b) reported that 2,4-D degradation rates in soils remained relatively constant with and without sunlight, suggesting that photodegradation is not an important process in the field. 2,4-D half-lives are short, ranging from a few days to several months but detectable residues can persist for up to a year (McCall et al., 1981). Since field studies to evaluate environmental fate and behavior are very expensive and provide a site-specific evaluation of environmental exposure, interest has developed in the use of models in predicting environmental fate for a range of circumstances (e.g., Jarvis et al., 2000; Zand-Parsa et al., 2006; Kandelous and Simunek, 2010; Bonfante et al., 2010; Tafteh and Sepaskhah, 2012; Moriasi et al., 2013). While certain problems may be solved using relatively simple analytical or semi-analytical models, other problems may require more sophisticated numerical models, either one- or multidimensional, that simulate water flow, solute transport, and a range of biogeochemical reactions (Simunek et al., 2008; Simunek, 2005; Batu, 2005). Acceptance by regulatory officials and the ability of the models to accurately describe environmental processes are two factors which should be considered in order to select the suitable model for simulating pesticide’s movement (Trevisan et al., 2000). PRZM has been validated in more than one study and it has a history of producing results acceptable to scientists and regulatory authorities (Cohen et al., 1995). PRZM3 (Pesticide Root Zone Model) is developed for the U.S. Environmental Protection Agency over a period of years (the first manual released in 1984). PRZM3 is the most recent version of a modeling system that links two subordinate models, PRZM and VADOFT, in order to predict pesticide transport and transformation down through the crop root and unsaturated zone. PRZM3 is a specific model for simulating pesticides as it is capable of simulating transport and transformation of the parent compound and as many as two daughter species. PRZM has been used in several studies by Young and Carleton (2006), Trevisan et al. (2000), Farenhorst et al. (2009), Noshadi et al. (2011), however, in some cases it is stated that the model does not accurately describe the behavior of applied pesticides (Banton and Villeneuve, 1989; Trevisan et al., 2000; Young and Carleton, 2006; Tahir and Lin, 2010). This can be due to complexity and diversity of soil hydraulic properties, which is inevitable. Besides, inaccuracy in predicting pesticide behavior can be explained by some simplifications assumed in the PRZM3 code. Since in unsaturated porous media water content impacts directly on solute transfer (Latrille, 2013), and the PRZM3 code uses the Richard’s equation (RE) for flow in the unsaturated zone, any simplification either in using Richard’s equation or in its numerical solution leads to inaccuracy in simulation. From the point of view of the PRZM’s developers, PRZM3 has two major deficiencies. First, in the PRZM subroutine Richard’s equation is simplified to simple drainage rules for simulating water movement in soil profile. This simplification leads to both inaccuracies in the final simulation results and increase in water balance error. When two component models (PRZM and VADOFT) link together, the second deficiency in

39

the model appears, which is due to the incompatibility between the hydraulics in the PRZM and VADOFT (Suarez, 2005). While VADOFT for movement of water in a saturated medium solves the Richards’ equation, PRZM simplifies Richards’ equation to the simple drainage rules to move water through the soil profile. Because of this incompatibility, there might be some times in which the value simulated for water content by PRZM is too high for VADOFT. The result would be water ponded at the interface which could not be attributed to PRZM or VADOFT (Suarez, 2005). According to Suarez (2005), this issue is very likely to happen in agricultural soils, where permeability of sub soils are typically lower than those of the root zone. Researchers solved RE analytically (Srivastava and Yeh, 1991; Lu and Zhang, 2004; Menziani et al., 2007; Tracy, 2007). These solutions are valid for simple conditions. Because unsaturated flow equation (RE) is highly nonlinear, an analytical solution is not possible except for special cases. Therefore, numerical methods such as finite difference or finite element were typically used to solve the unsaturated flow equation in the last thirty years (e.g., Haverkamp et al., 1977; Simunek et al., 1999; Rathfelder and Abriola, 1994). These solutions often suffer to some degree from mass balance errors as well as from numerical instability (Huang et al., 1994). In this research, we eliminated simplifications assumed in the original code of PRZM3 regarding RE and water movement equations and resolved the equations with an accurate numerical solution using MC-Cormack method which is numerically stable and avoids mass-balance errors. In addition, since the soil water has often been thought of as comprising two phases, one mobile and the other effectively immobile (Van Genuchten and Wierenga, 1976), we divided the soil water into mobile and immobile domains (MIM). We evaluated the prediction capability of the modified PRZM3 (in both modifications) using a data set describing results from field and laboratory experiments performed at the College of Agriculture of Shiraz University (Noshadi et al., 2011) and compared the results to those of the original PRZM3. We used inverse method in order to find unknown coefficients in MIM (mobileimmobile) modification including immobile soil water content and first rate coefficient. To calculate these coefficients Shuffled Complex Evaluation code, which is an optimization method, was added to the main program. 2. Materials and method 2.1. Water flow Darcy’s approach to groundwater flow considers flow as diffusion, and this makes it applicable to both saturated and unsaturated types of flow (water movement). Richard’s equation is obtained by application of the unsaturated Darcy’s flow law and mass conservation law (Richards, 1931). The moisture content based (h-based) form of Richards’ equation may be expressed as: ∂ ∂ ∂h = C(h) = ∂t ∂t ∂z



K(h)





∂h +1 ∂z

−S

(1)

where h is the water pressure head [L],  is the volumetric water content [L−3 L−3 ], t is time [T], z is the vertical dimension [L] (assumed positive upward), K(h) is the unsaturated hydraulic conductivity, C(h) is the retention water capacity [L−1 ] and S is the sink term [L−3 L−3 T−1 ]. In many practical engineering problems, fractures and fissures determine an environment that cannot be accurately homogenized by a single set of porous environmental parameters. Therefore mobile-immobile (MIM) models have been extended and applied to such problems in order to improve approximations. MIM model assumes that the flow is restricted to pores between the aggregates,

40

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

or inter-aggregates pores. These models assume that a matrix, which consists of immobile water pockets, can exchange, retain and store water, but it does not permit convective flow. This conceptualization leads to a two-region, dual-porosity and MIM type model of flow and transportation (Phillip, 1968; Van Genuchten and Wierenga, 1976). These models partition the liquid phase into mobile (flowing, inter-aggregate),  mo , and immobile (stagnant, intra-aggregate),  im , regions:  = mo + im

(2)

The MIM formulation that is used here for variably saturated flow employs the Richards Eq. (1) to describe flow in the macro-pores, and a mass balance equation to describe moisture dynamics in the ˇ unek ˚ et al., 2003): matrix as follows (Sim ∂mo ∂ = ∂t ∂z

  K



∂h +1 ∂z

− Sm − w (3)

∂mo = −Sim + w ∂t where Sm and Sim are sink terms for both regions, and  w is the transfer rate for water from inter to intra-aggregate pores. The mass transfer rate,  w for water between mobile and immobile regions, has been assumed to be proportional to the difference in effective water contents of the two regions as cited in other dualporosity studies (Phillip, 1968; Simunek et al., 2001); it uses the first-order rate equation: w =

∂im = ω(Sem − Seim ) ∂t

(4)

where  im is the matrix water content, ω is a first-order rate coefficient (T−1 ), and Sem and Seim are effective fluid saturations of the fracture and matrix regions, respectively.

In order to discretize the equation, it is written in the h-base form. In the predictor step in the MC-Cormack technique, the provisional value of h (pressure head) at time level n + 1 using FTFS (forward in time, forward in space) is calculated as follows: (Note that in all discretization, the top and the bottom indexes denote values for time and spatial discretization, respectively.) (n+1)∗

hi



=

0.5 × t n

C(h)i × z



− K(h)i+1



×



K(h)i

∂u ∂u +a =0 ∂t ∂x ∗

] = uni − a Predictor step → [un+1 i

(5) t n − uni ) (u x i+1

uni + [un+1 ] i 2

∗

−a

(6)

t ∗ ∗ ([un+1 ] − [un+1 ] ) i i−1 2x (7)

n

hi − hi+1 +1 z

n 

hi+1 − hi+2 +1 z

−

0.5 × t n

C(h)i

n

(w )i + hni

(8)

where i is spatial mesh, t is time interval between two nodes, and (n + 1)* illustrates the predicted term. It should be noted that in the MC-Cormack technique, the predictor is considered as t/2. Con∗ tinuing solution in the corrector step, the predicted value [hn+1 ] i is corrected using BTBS (backward in time, backward in space) as follows: (n+1)∗

MC-Cormack’s technique is a variant of the Lax-Wendroff approach but is much simpler in its application. Like the LaxWendroff method, the MC-Cormack method is also an explicit finite-difference technique that is second-order-accurate in terms of space and time. The technique was introduced in 1969 (MacCormack, 1969), and it became the most popular explicit finite-difference method for solving fluid flow problems in the following 15 years (Anderson, 1994). Moreover, results obtained with the MC-Cormack method are perfectly satisfactory for many fluid flow applications (Anderson, 1994). Application of the MC-Cormack method to equations’ procedures is made in two steps; the first involves a predictor step, and the second is a corrector step. In the predictor step, a “provisional” value of an interested variable (e.g., “U” in Eq. (5)) at time level n + 1 is estimated using FTFS (forward in time, forward in space), and ∗ ] is corrected using in the corrector step the predicted value [un+1 i BTBS (backward in time, backward in space). When applied to a linear, first order partial differential equation (1), which can model the transport of a substance, e.g., a pollutant, in a uniform fluid flow that is moving with velocity a, the predictor corrector method is as follows:

Corrector step → un+1 = i

2.3. Discretization of Richards’s equation using the MC-Cormack technique



2.2. MC-Cormack technique

Wave equation →

where i and n denote spatial and time discretization values, respectively, u(t, x) represents the concentration value of the pollutant at time t and at spatial position x, a is a fixed, nonzero constant, t is the time interval between two nodes, and (n + 1)* illustrates the predicted term. ∗ The term [un+1 ] is a temporary “predicted” value of U at the i time level n + 1. The corrector equation provides the final value of U at the time level n + 1 (un+1 ). Note that in the predictor equai tion, a forward difference is used for ∂u/∂x, while in the corrector equation, a backward difference is used. This differencing can be reversed and in some problems it is advantageous to do so.

hn+1 = 0.5 hni + hi i

 A=

n K(h)i−1



+

t n

C(h)i × z

×A−

n

hi−2 − hi−1 +1 z

n − K(h)i



t



n

n (w )i C(h)i

n 

(9)

hi − hi+1 +1 z

The algorithm of MC-Cormack technique can be seen as a flowchart in Fig. 2. 2.4. PRZM3 changes 2.4.1. General modifications PRZM3 has several versions, but unfortunately updated software has not been released in recent years, so it does not execute perfectly in new operating systems such as Windows Vista or 7 and in some cases it does not execute at all. The version of PRZM3 model system in current distribution is that built with a Lahey/Fujitsu compiler. In this study all PRZM3 codes are brought into a new user-friendly compiler environment named Visual Fortran. In order to make other modifications, PRZM3 codes are also able to compile in other compilers that users may be working with. Thus, a general update is done in PRZM3 codes making it compatible with new computer systems including Microsoft Windows Vista and 7. 2.4.2. Specific modifications As mentioned before, PRZM3 links two main subordinate models, PRZM and VADOFT. PRZM code accounts for the fate of pesticides and nitrogen in crop root zones and VADOFT is a code for simulating moisture movement and solute transportation in the vadose zone. Therefore, in order to add new codes to original codes, changes must be applied in all subroutines. In this study, two modifications were added to the original PRZM3 codes. In the first

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

Fig. 1. Flow chart of the shuffled complex evolution (SCE) method.

modification, Richards’ equation was solved with the MC-Cormack scheme in the VADOFT code and equations associated with water movement in the PRZM code were rewritten using the MC-Cormack scheme. All changes were applied to all subroutines and they were made compatible with original codes. In the second modification, soil water divides into two domains: mobile and immobile. Similarly Richards’ equation and all other equations associated with water movement were solved by the MC-Cormack scheme. 2.4.3. Shuffled complex evolution A shuffled complex evolution (SCE) is an evolutionary algorithm that performs local and global searches and its strategy combines the particular strengths of controlled random search (CRS) algorithms with the concept of competitive evolution (Holland, 1992) and the newly developed concept of complex shuffling. The strategy is presented and illustrated in Fig. 1. 2.5. Field data Data for analysis of the modified model were taken from Noshadi et al. (2011), from tests that were conducted at the College of Agriculture of Shiraz University located about 15 km north of Shiraz, at an altitude of 1810 m in south of Iran. Some major physical and chemical properties of the soil are described in Table 1. Herbicide 2,4-D, was sprayed on August 12th, 2007 on the surface of the soil at the rate of 3.5 kg a.i./ha in a corn field and then soil 2,4-D residue samples were measured periodically 6, 11, 21, 35 and 53 days after application, down to 1 m soil depth at 10 cm increments. According to the research by Noshadi et al. (2011), the experiment was done as a completely randomized design with three replications and plot sizes of 114 m2 (12 m × 9.5 m). The method

Fig. 2. MC-Cormack flowchart.

41

42

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

Table 1 Characteristic of soil in the experimental site. Depth (cm)

Texture

Sand (%)

Silt (%)

Clay (%)

FCa (%)

PWPa (%)

OM (g kg−1 )

pHa

CECa (cmol kg−1 )

ECa (dS m−1 )

Ks (×10−6 cm s−1 )

b (g cm−3 )

0–10 10–20 20–30 30–40 40–50 50–60 60–70 70–80 80–90 90–100 100–110

SiC SiC SiC SiC SiC SiC SiCL SiCL SiCL SiCL SiCL

12.00 12.83 13.67 12.00 12.33 12.67 17.33 17.33 17.33 19.33 19.33

46.00 44.83 43.67 44.00 42.33 40.67 43.33 46.33 49.33 43.33 45.33

42.00 42.33 42.67 44.00 45.33 46.67 39.33 36.33 33.33 37.33 35.33

23 23 23 24 24 24 24 24 24 24 –

8 8 8 11 11 11 11 11 11 11 –

14.5 14.6 14.6 11.7 10.2 8.7 6.9 8.5 10.1 5.3 4.6

7.60 7.73 7.73 7.0 7.90 7.73 7.43 7.40 – – –

20.8 21.3 21.8 20.8 20.82 20.84 20.84 20.84 – – –

0.762 0.569 0.728 0.507 0.408 0.513 0.450 0.450 – – –

203 203 47.1 47.1 20.9 20.9 10.2 10.2 4.88 4.88 4.88

1.26 1.43 1.43 1.43 1.43 1.43 1.43 1.43 1.43 1.43 –

All parameters are an average of three replications. a FC, field capacity; PWP, permanent wilting point; OM, organic matter; Ks , saturated hydraulic conductivity; b , soil apparent density; SiC, silty clay; SiCL, silty clay loam. Table 2 Irrigation depth for normal and deficit irrigations during growth season. Date

Normal irrigation

7 Jun. 15 Jun. 22 Jun. 6 Jul. 14 Jul. 20 Jul. 31 Jul. 11 Aug. 30 Aug. 7 Sep. 23 Sep.

71 98.33 72.6 66.7 71.3 54.33 94.33 153.33 106.66 96.66 105

Sum

990

Deficit irrigation 64 88 54.4 46.41 40.66 54 63.33 96.66 84.33 70 73.33

Table 3 PRZM-3 main input parameters for simulating 2,4-D movement. Meteorological input parameters

Crop input parameters

735

of irrigation was that of a solid set sprinkler with two irrigation treatments of normal and deficit irrigation with irrigation depths of 910 and 735 mm, respectively (Table 2). Deficit irrigation was applied as reductions in irrigation time. Soil temperature was also recorded at the following depths: 5, 10, 30, 50 and 100 cm. Based on soil conditions such as temperature, moisture content and texture at different depths, herbicide behavior differs in the soil profile (Mersie and Foy, 1985). Therefore, soil depth was divided into 10 cm increments and soil samples were taken periodically with a handheld soil auger down to a 1-m depth at 10 cm increment from each plot and then transported to the laboratory and frozen at −20 ◦ C. 2.6. Modeling parameters Tables 3 and 4 show the main input parameters used in simulating 2,4-D movement in the soil profile. The input parameters listed in Tables 3 and 4 were the same in all modeling runs. It can be seen that, given exactly the same experimental data, different values can be chosen for some of the input parameters. For example, in Table 3, evapotranspiration value which is very important for the estimation of soil water content, can be calculated using monthly daylight data or a factor combined with the evapotranspiration data. We used daily pan evaporation data for this quantity. Volatilization of 2,4-D is expected to be very low based on its extremely low Henry’s Law constant (1.02 × 10−8 atm-cu m/mole) (USEPA, 1987). Diffusion coefficient and vaporization enthalpy of 2,4-D were assumed to be negligible and they were set to zero (USEPA, 1987). The number and the thickness of horizons and the thickness of horizon compartments could be entered differently. The choice of horizon thickness was made based upon the soil data provided with the data set. Considering soil hydraulic properties for each compartment, field capacity, permanent wilting point, saturated hydraulic conductivity, organic matter and soil apparent density were included in the input file. The values for the mentioned properties were

Management input parameters

Numerical discretization of the soil profile Hydraulic properties of soil profile

Pan factor Snowmelt factor Minimum depth of evaporation (cm) Pan factor flag (zero = daily pan evaporation data read; one = temperature data read)

0.75 0.00 10 0

Number of different crop Emergence date Maturation date Harvest date Maximum interception storage (cm) Maximum rooting depth (cm) Maximum areal coverage of the canopy (%) Surface condition after harvest residue Maximum canopy height at maturation (cm)

1 5/28/2007 11/10/2007 20/10/2007 30 90 80

Application date Depth of application (cm) Total application of 2,4-D (kg/ha) Total depth of soil core (cm) Diffusion coefficient in air (cm2 /day) Henry’s law constant (atm-cu m/mole) Enthalpy of vaporization (kcal/mol)

12/7/2007 4 3.5 90 0 1.02 × 10−8

Total number of horizons Horizon thickness (cm) Compartment thickness (cm)

9 10 2

Field capacity Permanent wilting point Saturated hydraulic conductivity Organic matter Soil apparent density

Listed in Table 1 for each soil layer

Residue 200

0

Table 4 Decay rate (per day) in dissolved and adsorbed phase and partition coefficient (KD) (cm−3 g−1 ) of 2,4-D for each horizon used in calibrations and simulations. Horizon

Decay rate (per day) in dissolved phase

Decay rate (per day) in adsorbed phase

Partition coefficient (cm−3 g−1 )

1 2 3 4 5 6 7 8 9

0.098 0.058 0.049 0.021 0.021 0.021 0.021 0.021 0.021

0.098 0.058 0.049 0.021 0.021 0.021 0.021 0.021 0.021

1.61 1.60 1.60 1.60 1.60 1.60 1.60 1.60 1.60

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

43

Table 5 Statistical parameters in normal irrigation treatment. Days after application

NRMSE

CRM

d

PRZM3

PRZM3-MC

PRZM3-MC-MIM

PRZM3

PRZM3-MC

PRZM3-MC-MIM

PRZM3

PRZM3-MC

PRZM3-MC-MIM

8 13 23 30 37 57

0.39 0.13 0.26 0.59 1.03 1.11

0.14 0.16 0.14 0.22 0.40 0.60

0.15 0.13 0.13 0.24 0.33 0.39

−0.36 −0.09 0.18 0.44 0.42 0.20

−0.09 −0.03 −0.09 −0.06 −0.15 0.17

−0.12 0.00 −0.05 −0.01 −0.18 −0.01

0.95 0.96 0.49 0.30 0.49 0.78

0.99 0.95 0.33 0.61 0.90 0.96

0.99 0.97 0.72 0.64 0.89 0.94

Mean

0.59

0.28

0.23

0.13

−0.04

−0.06

0.66

0.79

0.86

measured and listed in Table 1 for each compartment. Table 4 shows the decay rates and sorption coefficients. PRZM-3 requires only the input of linear sorption coefficients and a first-order decay rate. The value of these parameters which are vital to the prediction of leachate concentration were selected during calibration process. To calibrate a parameter an initial value for the given soil type (based on the range of the parameter found in the literature) was assumed in the input file. After several runs, the initial value was changed in its logical range in order to find the best fit between the simulated and measured data. When the minimum difference between simulated and measured data is attained the parameter is calibrated. It should be noted that, two parameters including the immobile water content and ω (the first-order rate coefficient) should be given as program input when users run the model with the MIM modification. The model uses the SCE algorithm and calculates optimized values of the mentioned parameters. 2.7. Statistical analysis For testing 2,4-D simulation output, the following statistical parameters were used to make comparisons between simulated data using original PRZM and modified PRZM (MC-Cormack and MIM model) with data measured from the field (Noshadi et al., 2011):

  n 100  1 NRMSE = (Pi − Oi )2 ¯ O

n CRM =

d=1−

n

Oi −

n

n

i=1

(10)

i=1

P i=1 i

(11)

O i=1 i

n (Pi − Oi )2

n i=1 ∗ ∗ 2

¯ Pi∗ = Pi − O

i=1

(|Pi | − |Oi |)

(12)

¯ Oi∗ = Oi − O

   Valuesimulated − ValueMeasured   × 100 Value

% Error = 

(13)

Measured

In the above equations, Oi and Pi are observed and predicted 2,4-D concentrations, respectively, n is the number of observations and ¯ represents the mean of observed data values. O NRMSE (normal root mean square error) provides the total difference between measured and simulated data proportioned against means of measured data values. The lower limit for NRMSE is zero, which occurs when there is no difference between such paired data. Obviously, a small value of NRMSE indicates a more accurate simulation. CRM (coefficient of residual mass) is an indication of consistent error in distribution of all simulated values across all measurements with no consideration of the order of measurements. A value of zero for CRM indicates no bias in distribution of simulated values with respect to measured values. The index of agreement, d, was calculated for assessing accuracy of simulated data. The maximum value for d is one, which occurs when simulated values are completely identical to measured values. Percentage error is defined as a percentage of the difference between an approximate or simulated value and an exact or measured value. A lower value for percentage error means that the simulation result is closer to the measured data.

3. Results and discussion In this study, simulation results of the modified PRZM3 (8, 13, 30 and 57 days after 2,4-D application) were compared to results of Noshadi et al. (2011) in which 2,4-D concentrations in the soil profile were measured under two irrigation treatments (normal irrigation and deficit irrigation). Values of d, NRMSE and CRM are shown in Tables 5 and 6. The mean of these parameters for PRZM3 in normal irrigation for all measured and simulated data were 0.66, 0.59 and 0.13, for simulated data using MC-Cormack (PRZMMC) were 0.79, 0.28 and −0.04 and for mobile-immobile approach (PRZM-MC-MIM) were 0.86, 0.23 and −0.06, respectively. In deficit irrigation d, NRMSE and CRM for all measured and simulated data were 0.65, 0.52 and 0.08, for simulated data using the MC-Cormack technique the values were 0.77, 0.38 and −0.24, and for mobileimmobile approach they were 0.73, 0.36 and −0.24, respectively. In addition, to evaluate the accuracy of modified models, percentage error can be used through Eq. (13). The average of percentage error obtained for all simulation results in the original

Table 6 Statistical parameters in deficit irrigation treatment. Days after application

NRMSE

CRM

d

PRZM3

PRZM3-MC

PRZM3-MC-MIM

PRZM3

PRZM3-MC

PRZM3-MC-MIM

PRZM3

PRZM3-MC

PRZM3-MC-MIM

8 13 23 30 37 57

0.48 0.30 0.40 0.55 0.54 0.86

0.18 0.20 0.41 0.60 0.54 0.38

0.21 0.18 0.30 0.43 0.69 0.38

0.46 0.28 0.02 0.24 0.30 0.70

−0.18 −0.19 −0.36 −0.52 −0.41 0.21

−0.21 −0.15 −0.29 −0.32 −0.56 0.10

0.91 0.90 0.71 0.49 0.49 0.42

0.99 0.96 0.54 0.50 0.78 0.87

0.98 0.97 0.69 0.60 0.42 0.74

Mean

0.52

0.38

0.36

0.08

−0.24

−0.24

0.65

0.77

0.73

44

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

Fig. 3. Measured and PRZM3 simulation of 2,4-D in different soil depths for normal irrigation: (a) 8 days after application, (b) 13 days after application, (c) 30 days after application and (d) 57 days after application.

model was 44.88% while this reduced to 27.05% using MC-Cormack and 23.62% in PRZM-MC-MIM. A further point to consider is that in the original PRZM3, water balance error increased noticeably, particularly with an increase in the number of simulation days. In this research, the cumulative amount of water balance error reached 35 × 10−3 in the 56th day. This amount decreased remarkably following model modification, being recorded at 21 × 10−5 using MC-Cormack and 44 × 10−5 in PRZM-MC-MIM. Following the MIM modification, SCE code was used by the modified model to optimize the MIM unknown coefficients. It was found that the best fit between the simulated and measured data was attained when immobile soil water content was 0.08 cm3 cm−3 and the first rate coefficient was 0.269 (min−1 ). Note that, the initial value of immobile soil water content was set to 0.1 and the first rate coefficient was set to 0.5. 3.1. Normal irrigation On July 20th, 8 days after 2,4-D application, the maximum depth at which 2,4-D was detected, was 30 cm (Noshadi et al., 2011). The original model was able to simulate 2,4-D concentration in the first and second layers relatively accurately, although the original model’s simulation for the third layer was far different from that of the field data. For this date, d, NRMSE and CRM were 0.94, 0.39, −0.36, respectively, which shows relatively good agreement between measured and simulated data. For this date after modification of the model, using the MC-Cormack approach (PRZM3-MC), results shifted to better alignment with the field data, especially in the third layer

which simulated 2.66 mg kg−1 soil, compared to the original model, which simulated 5.02 mg kg−1 soil. A reduction in NRMSE to 0.13 and an increase of d value to 0.99 confirms this claim. This can also be seen in the MIM approach (PRZM3-MC-MIM), although 2,4D concentration was slightly higher for the three layers involved. Values of d, NRMSE and CRM were 0.99, 0.15 and −0.12, respectively, which demonstrates the best agreement between measured and simulated data for this date among models (Fig. 3a). Average concentrations of 2,4-D in 0–10, 10–20 and 20–30 cm soil layers were reduced to 8, 3.5 and 4.86 mg kg−1 soil, on August 3rd, due to degradation 13 days after adding 2,4-D. Values of d, NRMSE and CRM were 0.96, 0.13 and −0.09, respectively showing that simulations were accurate for the original PRZM3. For this date, after modification and using the MC-Cormack approach, results for the first two layers did not show significant improvement, although data simulation in the third layer was in better agreement with that of the field data. In the MIM approach, a notable difference was determined compared to MC-Cormack approach, in the second layer where 2,4-D concentration was simulated to be 3.62 mg kg−1 soil. Based on the calculations for the MC-Cormack approach, d, NRMSE and CRM were 0.95, 0.16 and −0.03, and for the MIM approach these were 0.97, 0.13 and 0.003, respectively, showing better simulation in comparison to the original model (Fig. 3b). Thirty days after application, 2,4-D leached to the depth of 40 cm. Concentrations reduced in the first and second layers; however, higher concentrations were recorded in the fourth layer. This may be due to an application of 153.33 mm irrigation water on August 11, which was higher than other irrigation applications due to more demand for water by the corn crop. Predictions for the first three layers were relatively acceptable. In the fourth layer, the

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

45

Fig. 4. Measured and PRZM3 simulation of 2,4-D in different soil depths for deficit irrigation: (a) 8 days after application, (b) 13 days after application, (c) 30 days after application and (d) 57 days after application.

measured value was 3.27, whereas 0.33 mg kg−1 soil was simulated for this layer with the original model. With the exception of the third layer, simulation after modification illustrated a considerably better agreement. This agreement was more significant in the fourth layer where 2,4-D concentration were calculated to be 2.52 and 2.63 mg kg−1 soil with MC-Cormack and MIM approaches, respectively. Values of d, NRMSE and CRM for original PRZM were 0.3, 0.59 and 0.44, respectively, while these statistical parameters for modified models indicated general improvements of simulations. For the MC-Cormack approach, d, NRMSE and CRM were 0.61, 0.22 and −0.06, and for MIM approach they were 0.64, 0.24 and −0.0, respectively (Fig. 3c). Finally, fifty-seven days after application, 2,4-D concentrations in the first five 10 cm layers of the soil were 0, 0, 0, 3.2 and 0 mg kg−1 soil. The model was only able to make good predictions of 2,4-D concentrations in the fourth layer; however, over-estimation occurred in the other layers. Values of d, NRMSE and CRM were 0.78, 1.11 and 0.2, respectively, showing that simulation was not satisfactory. Over-estimation in the first three layers happened in the modified model as well, but it was less evident. Simulated 2,4-D concentration for third layer was 0.23 and 0.37 mg kg−1 soil for the MC-Cormack and MIM approaches, respectively, while it was predicted to be 0.76 mg kg−1 soil by original PRZM3. Besides, 2.27 and 2.71 mg kg−1 soil of 2,4-D concentration in the fourth layer, calculated by the MC-Cormack and MIM approaches, showed good agreement with measured data. Based on these calculations, for the MC-Cormack approach, d, NRMSE and CRM values were 0.96, 0.6 and 0.17 respectively, and for MIM approach these values were 0.94, 0.39 and −0.01, respectively (Fig. 3d).

The average percentage error obtained for all simulation results in the original model was 43.94%, while this reduced to 22.82% using MC-Cormack and 20.08% in PRZM-MC-MIM. 3.2. Deficit irrigation As it can be seen in Fig. 4a, 8 days after adding 2,4-D, the pesticide reached the third layer (Noshadi et al., 2011). Note that less soil moisture causes less microbial activity and therefore 2,4-D degradation decreased (Cattaneo et al., 1997). The original model could not accurately simulate 2,4-D concentrations for this date, especially the third layer in which simulated concentration was three times higher than that of the field data. For the original model, d, NRMSE and CRM were 0.91, 0.48 and −0.46, respectively. For this date after modification using the MC-Cormack modification (PRZM3-MC), better results were achieved, especially in the third layer which 2,4-D concentration was simulated to be 3.23 mg kg−1 soil, compared to the original model, which was simulated to be 6.08 mg kg−1 soil. A reduction in NRMSE to 0.18 and an increase of d value to 0.98 confirm this claim. Improved results can also be seen in the MIM approach, although 2,4-D concentration was slightly higher for the three layers involved. For the MIM approach, values of d, NRMSE and CRM were 0.98, 0.21 and −0.21, respectively. As Fig. 4b demonstrates, average concentrations of 2,4-D in 0–10, 10–20 and 20–30 cm soil layers were reduced to 9.5, 4 and 3.2 mg kg−1 soil, on August 3rd. Values of d, NRMSE and CRM were 0.9, 0.3 and −0.28, respectively showing that the simulations using the original model were not accurate enough. For this date after modification and using the MC-Cormack approach, results for the first two layers did not show significant improvement, although in

46

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47

the third layer simulation data were in better agreement with field data. In the MIM approach, a notable difference was determined in the second layer where 2,4-D concentration was simulated to be 4.35 mg kg−1 soil compared to the MC-Cormack approach. Based on the calculations for the MC-Cormack approach, d, NRMSE and CRM were 0.96, 0.2 and −0.19, and for the MIM approach these parameters were 0.97, 0.18 and −0.15, respectively, showing better simulation in comparison to the original model (Fig. 4b). Thirty days after application, 2,4-D was recorded in the depth of 40 cm. On this date, simulations for the first three layers were relatively acceptable. In the fourth layer, the measured value was 3.1, whereas 0.45 mg kg−1 soil was simulated with the original model. Unlike other simulations in this date, with the exception of the fourth layer, simulations after modification with the MCCormack approach did not illustrate better agreement; however, in the MIM approach simulations shifted to better results. The values of d, NRMSE and CRM for the original model were 0.49, 0.55 and 0.24, respectively, while for the MC-Cormack approach the statistical parameters were 0.50, 0.60 and −0.52, respectively, and for the MIM approach were 0.60, 0.43 and −0.32, respectively (Fig. 4c). The last measurements were carried out on September 14th (fifty-seven days after application) when 2,4-D concentrations in the first four layers of the soil were 1.8, 1.43, 1.1 and 3.3 mg kg−1 soil. 2,4-D concentrations simulated with the original model in the layers were less than observed values. Values of d, NRMSE and CRM were 0.42, 0.86 and 0.7, respectively, showing that simulations with the original model were not satisfactory. In the MC-Cormack approach, however, simulated concentrations were closer to the observed data especially in the first and the last layers where 2,4-D concentrations were simulated to be 2 and 3.07 mg kg−1 soil, respectively. In the MIM approach, simulation results were also satisfactory. Based on the calculations for the MC-Cormack approach, d, NRMSE and CRM were 0.87, 0.38 and 0.21, and for the MIM approach these were 0.74, 0.38 and 0.1, respectively (Fig. 4d). The average percentage error obtained for all simulation results in the original model was 53.81%, while this reduced to 31.28% using MC-Cormack and 27.16% in PRZM-MC-MIM.

4. Conclusion In this paper PRZM3 was modified by supplementing two new modifications: the MC-Cormack technique and mobileimmobile (MIM) in which Richards’ equation was solved using the MC-Cormack technique and was changed to the form of a mobile-immobile evaluation. All equations related to moisture and RE were changed in the code. In addition, shuffled complex evolution was added to codes in order to calculate mobile-immobile coefficients. Statistical analysis including NRMSE, CRS, d and percentage error demonstrated significant improvements in simulation results after modification with both MC-Cormack and MIM approaches. In addition, water balance error decreased noticeably after modification which confirmed the improvement of the model. Although the Mac-Cormack technique used in this study was found to be numerically stable and led to better simulations, model running time increased 1.5 times. This was due to using two time steps in equation solving using MC-Cormack method and addition of shuffled complex method to the model codes. Moreover, input parameters increased in number due to modifying the original model to mobile-immobile approach. In the PRZM-MC-MIM modified model, users need to enter immobile water content and ω (first-order rate coefficient) values in addition to the input parameters in the original model. As measuring these parameters is

complicated, care should be taken prior to running PRZM-MC-MIM model (users, however, can apply inverse method using shuffled complex evaluation to calculate the values of such parameters). References Anderson Jr., J.D., 1994. Computational Fluid Dynamics: The Basics with Applications, 1st ed. McGraw Hill, New York. Banton, O., Villeneuve, J.P., 1989. Evaluation of groundwater vulnerability to pesticides: a comparison between the pesticide drastic index and the PRZM leaching quantities. J. Contam. Hydrol. 4 (3), 285–286. Batu, V., 2005. Applied Flow and Solute Transport Modeling in Aquifers: Fundamental Principles and Analytical and Numerical Methods. CRC Press, Boca Raton, FL, pp. 720 pp. Bonfante, A., Basile, A., Acutis, M., De Mascellis, R., Manna, P., Perego, A., Terribile, F., 2010. SWAP, CropSyst and MACRO comparison in two contrasting soils cropped with maize in Northern Italy. Agric. Water Manage. 97, 1051–1062. Cattaneo, M.V., Masson, C., Greer, W.C., 1997. The influence of moisture on microbial transport, survival and 2,4-D biodegradation with a genetically marked Burkholderia capacious in unsaturated soil. Biodegradation 8, 87–96. Cohen, S.Z., Wauchope, R.D., Klein, R.D., Eadsforth, C.V., Graney, R., 1995. Offsite transport of pesticides in water: mathematical models of pesticide leaching and runoff. J. Pure Appl. Chem. 67, 2109–2148. Farenhorst, A., McQueen, D.A.R., Saiyed, I., Hilderbrand, C., Li, S., Lobb, D.A., Messing, P., Schumacher, T.E., Papiernik, S.K., Lindstrom, M.J., 2009. Variations in soil properties and herbicide sorption coefficients with depth in relation to PRZM (pesticide root zone model) calculations. Geoderma 150 (3–4), 267–277. Foster, R.K., McKercher, R.B., 1973. Laboratory incubation studies of chlorophenoxyacetic acids in chernozemic soils. Soil Biol. Biochem. 5, 333–337. Halter, M., 1980. 2,4-D in the aquatic environment. Section II. In: Shearer, R., Halter, M. (Eds.), Literature Reviews of Four Selected Herbicides: 2,4-D, Dichlobenil, Diquat & Endothall. Study conducted for Seattle, Washington METRO. Han, S.O., New, P.B., 1994. Effect of water availability on degradation of 2,4dichlorphenoxyacetic acid (2,4-D) by soil microorganisms. Soil Biol. Biochem. 26 (12), 1689–1697. Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J., Vachaud, G., 1977. A comparison of numerical simulation models for one dimensional infiltration. Soil Sci. Soc. Am. J. 41 (2), 285–294. Helling, C.S., Kearney, P.C., Alexander, M., 1971. Behavior of pesticides in soil. Adv. Agron. 23, 147–240. Hemmett, R.B., Faust, S.D., 1969. Biodegradation kinetics of 2,4dichlorophenoxyacetic acid by aquatic microorganisms. Residue Rev. 29, 191–207. Holland, J.H., 1992. Adaptation in Natural and Artificial Systems, 1st ed. A Bradford Book, London. Huang, K., Zhang, R., Van Genuchten, M.T., 1994. An Eulerian–Lagrangian approach with an adaptively corrected method of characteristics to simulate variably saturated water flow. Water Resour. Res. 30 (2), 499–507. Jarvis, N.J., Brown, C.D., Granitza, E., 2000. Sources of error in model predictions of pesticide leaching: a case study using the MACRO model. Agric. Water Manage. 44, 247–262. Johnson, W.G., Lavy, T.L., Gbur, E.E., 1995a. Sorption, mobility, and degradation of triclopyr and 2,4-D and four soils. Weed Sci. 43, 678–684. Johnson, W.G., Lavy, T.L., Gbur, E.E., 1995b. Persistence of triclopyr and 2,4-D in flooded and nonflooded soils. J. Environ. Qual. 24, 493–497. Kandelous, M.M., Simunek, J., 2010. Numerical simulation of water movement in a subsurface drip irrigation system under field and laboratory condition using HYDRUS-2D. Agric. Water Manage. 97 (7), 1070–1076. Latrille, C., 2013. Effect of water content on dispersion of transferred solute in unsaturated porous media. Procedia Earth Planet. Sci. 7, 463–466. Lu, Z., Zhang, D., 2004. Analytical solution to steady state unsaturated flow in layered randomly heterogeneous soils via Kirchhoff transformation. Adv. Water Resour. 27 (8), 775–784. MacCormack, R.W., 1969. The effect of viscosity in hypervelocity impact cratering. AIAA Paper 69-354. McCall, P.J., Vrona, S.A., Kelley, S.S., 1981. Fate of uniformly carbon-14 ring labeled 2,4,5-trichlorophenoxyacetic acid and 2,4-dichlorophenoxyacetic acid. J. Agric. Food Chem. 29, 100–107. Menziani, M., Pugnaghi, S., Vincenzi, S., 2007. Analytical solutions of the linearized Richards equation for discrete arbitrary initial and boundary conditions. J. Hydrol. 332, 214–225. Mersie, W., Foy, C.L., 1985. Phytotoxicity and adsorption of chlorsulfuron as affected by soil properties. Weed Sci. 33, 564–568. Moriasi, D.N., Gowda, P., Arnold, J.G., Mulla, D.J., Ale, S., Steiner, J.L., 2013. Modeling the impact of nitrogen fertilizer application and tile drain configuration on nitrate leaching using SWAT. Agric. Water Manage. 130, 36–43. Noshadi, M., Foroharfar, F., Amin, S., Mafton, M., 2011. Measuring and simulating 2,4-D residues in silty clay soil profile under two water regimes using LEACHP model. Iran. Agric. Res. 30, 33–46. Oh, K.-H., Tuovinen, O.H., 1991. Bacterial degradation of phenoxy herbicide mixtures 2,4-D and MCPP. Bull. Environ. Contam. Toxicol. 47, 222–229. Ou, L.-T., 1984. 2,4-D degradation and 2,4-D degrading microorganisms in soils. Soil Sci. 137 (2), 100–107.

M. Noshadi, S. Jamshidi / Agricultural Water Management 143 (2014) 38–47 Phillip, J.R., 1968. The theory of absorption in aggregated media. Aust. J. Soil Res. 6, 1–19. Que Hee, S.S., Sutherland, R.G., 1981. The Phenoxyalkanoic Herbicides, vol. I: Chemistry, Analysis, and Environmental Pollution. CRC Press, Inc., Boca Raton, Florida, pp. 319 pp. Rathfelder, K., Abriola, L.M., 1994. Mass conservative numerical solutions of the head based Richards equation. Water Resour. Res. 30 (9), 2579–2586. Richards, L.A., 1931. Capillary conduction of liquids through porous medium. Physics 1, 318–333. Sandmann, E.R.I.C., Loos, M.A., van Dyk, L.P., 1988. The microbial degradation of 2,4-dichlorophenoxyacetic acid in soil. Rev. Environ. Contam. Toxicol. 101, 1–53. Shaw, L.J., Burns, R.G., 1998. Biodegradation of 2,4-D in a noncontaminated grassland soil profile. J. Environ. Qual. 27, 1464–1471. Simunek, J., 2005. Models of water flow and solute transport in the unsaturated zone. In: Encyclopedia of Hydrological Sciences, 1st ed. Wiley, California. ˇ Simunek, J., Van Genuchten, M.Th., Sejna, M., 2008. The HYDRUS-1D Software Package for Simulating the One-dimensional Movement of Water, Heat, and Multiple Solutes in Variably-saturated Media. Version 4.0. HYDRUS Software Ser. 3. Dep. of Environmental Sciences, Univ. of California, Riverside. Simunek, J., Wendroth, O., Van Genuchten, M.Th., 1999. Soil hydraulic properties from laboratory evaporation experiments by parameter estimation. In: Van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Proceedings of the International Workshop, Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. University of California, Riverside, CA, pp. 713–724. Simunek, J., Wendroth, O., Wypler, N., van Genuchten, M.Th., 2001. Nonequilibrium water flow characterized from an upward infiltration experiment. Eur. J. Soil Sci. 52 (1), 13–24. ˇ unek, ˚ Sim J., Jarvis, N.J., van Genuchten, Th.M., Gärdenäs, A., 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272, 14–35. Smith, A.E., Aubin, A.J., 1994. Loss of enhanced biodegradation of 2,4-D and MCPA in a field soil following cessation of repeated herbicide applications. Bull. Environ. Contam. Toxicol. 53, 7–11.

47

Srivastava, R., Yeh, T.C.J., 1991. Analytical solutions for one dimensional, transient infiltration toward the water table in homogeneous and layered soils. Water Resour. Res. 27 (5), 753–762. Suarez, A.L., 2005. PRZM-3, A Model for Predicting Pesticide and Nitrogen Fate in the Crop Root and Unsaturated Soil Zones: User’s Manual for Release 3.12.2. National Exposure Research Laboratory, U.S. Environmental Protection Agency, Athens, GA 30605-2700. Tafteh, A., Sepaskhah, A.R., 2012. Application of HYDRUS-1D model for simulating water and nitrate leaching from continuous and alternate furrow irrigated rapeseed and maize fields. Agric. Water Manage. 113, 19–29. Tahir, A., Lin, A.H., 2010. GIS based ArcPRZM-3 model for bentazon leaching towards groundwater. J. Environ. Sci. 22 (12), 1854–1859. Tracy, T.F., 2007. Three-dimensional analytical solutions of Richards’ equation for a box-shaped soil sample with piecewise-constant head boundary conditions on the top. J. Hydrol. 336, 391–400. Trevisan, M., Errera, G., Goerlitz, B., Remy, B., Sweeney, P., 2000. Modelling ethoprophos and bentazone fate in a sandy humic soil with primary pesticide fate model PRZM-2. Agric. Water Manage. 44, 317–335. USEPA, 1987. Office of Ground Water and Drinking Water. National Primary Drinking Water Regulations, Technical Fact Sheet on 2,4-D. Available: http://www.epa.gov/ogwdw/pdfs/factsheets/soc/tech/24-d.pdf Van Genuchten, M., Wierenga, P.J.Th., 1976. Mass transfer studies in sorbing porous media: I. Analytical solutions. Soil Sci. Soc. Am. J. 40, 473–481. Wang, Y.-S., Yen, J.-H., Hsieh, Y.-N., Chen, Y.-L., 1994. Dissipation of 2,4-D, glyphosate, and paraquat in river water. Water Air Soil Pollut. 72, 1–7. Wilson, R.D., Geronimo, J., Armbruster, J.A., 1997. 2,4-D dissipation in field soils after applications of 2,4-D dimethylamine salt and 2,4-D 2-ethylhexyl ester. Environ. Toxicol. Chem. 16 (6), 1239–1246. Worthing, C.R., Hance, R.J., 1991. The Pesticide Manual. BCPC, Surrey, UK. Young, D.F., Carleton, J.N., 2006. Implementation of a probabilistic curve number method in the PRZM runoff model. Environ. Model. Softw. 21, 1172–1179. Zand-Parsa, Sh., Sepaskhah, A.R., Ronaghi, A., 2006. Development and evaluation of integrated water and nitrogen model for maize. Agric. Water Manage. 81, 227–256.