Numerical Simulation of Pesticide Transport and Fate for Water Management in the Fucino Plain, Italy ()
1. Introduction
Detection of pesticides and their breakdown products in surface water and in groundwater systems is needed due to their detrimental environmental impact. Pesticides can reach surface water by runoff or by infiltration and seepage into water courses. Infiltration of rainfall and irrigation waters can leach pesticides to the water table, causing contamination of aquifers. Prevalence of one of these processes mainly depends on the amount and distribution of recharge and on the hydraulic conductivity of the unsaturated zone.
The transport and fate of contaminants is a process that requires the knowledge of many parameters, such as soil physical, (bio)chemical and hydrogeological properties as well as climatic conditions and agricultural practices. The types of pesticide (herbicide, insecticide, fungicide) and crop also influence the process. Due to the high number of variables, a study of contamination by pesticides is not feasible by large-scale in situ experiments. A different approach is the study of entire catchment based on observation points [1–4], taking into account both hydrogeological context and hydrodynamics [5–7]. An easier way is represented by laboratory experiments [8–12], while a credible and diffused tool is the use of mathematical models to simulate the transport and fate process. Among the large number of transport models available [13–26], the integrated pesticide transport model (IPTM) developed by [27] is chosen for this study. Thus, this paper is an application of a well-tested and known pesticide transport model [27,28].
To compare simulation results with field data, use was made of field data collected from the agricultural area of the Fucino Plain, Central Italy, which has been studied in the last ten years from hydrological, hydrogeological, and agricultural points of view [29–32].
The Fucino Plain (200 km2 wide) was the largest lake in Central Italy (Figure 1), reclaimed in the 1800s for agricultural purposes. The fractured carbonate aquifers surrounding the Plain feed high-discharge springs and streambed springs, which ensure steady discharges even during the dry season [30]. The heterogeneous aquifer of the Plain, having a variable vertical permeability, is supplied by groundwater seepage and by direct infiltration from rainfall. The long-term water balance of the Plain has 700 mm/y of rainfall, 450 mm/y of evapotranspiration,
Figure 1. Hydrological setting of the Fucino Plain. Main network of artificial canals is shown by black thin lines. 1: Plain aquifer, corresponding to the agricultural area; 2: fan and detrital deposits connecting carbonate aquifers to the Plain aquifer; 3: ancient alluvial deposits (aquitard); 4: terrigenous deposits (regional aquiclude); 5: carbonate aquifers (recharge area of springs); 6: climatic stations; 7: main springs; 8: main streambed springs; 9: public irrigation wells; 10: drinking water wells; 11: main directions of groundwater flow.
and consequently 250 mm/y of water excess in the October-March period [30]. As a consequence, crop irrigation was based on a sustainable use of surface water and, increasingly, groundwater through the 1980s. However, during the past 15 years, due to the concurrence of a natural decrease in spring discharge and the increase of pumping and water requirements for irrigation, a significant water shortage in the summer has been observed. Ignoring signs of water and environmental imbalance of the water-man-agriculture system, historical farm crops, mostly wheat, potato and sugar beet, were progressively and rapidly replaced by water-intensive horticultural crops [30]. This transition to an intensive agriculture has been accompanied by high water demand and by wide use of pesticides. In addition, the decreasing discharge in the canals nullifies the effect of dilution of the pollutants entering surface waters, mainly residues of fertilizers [32] and pesticides. In this framework, surveys conducted on both surface water and groundwater have shown the presence of pesticides, mainly Linuron, Dichloran, and Carbaryl, with concentrations ranging from 0.02 to 2.8 g/L in surface water and from 0.03 to 0.5 g/L in groundwater [33]. Those concentrations are frequently higher than threshold values (0.1 g/L in groundwater and surface water for human use) allowed by Italian laws and European Union Directives [34–36].
There is a need to clarify how surface runoff and/or infiltration processes allow pesticides to reach canals and/or groundwater. Because of focus on pesticide occurrence in surface waters and groundwater, gaseous and adsorbed components of pesticide fate are not studied in detail. The IPTM-CS model [28,37] takes into account those components as functions of parameters that have been considered constant during these simulations (e.g., organic matter content for the adsorbed component). Contribution of pesticides on volatilization, adsorption, and decay has been considered as losses for water systems.
The aim of this study is to verify by numerical simulation modeling how various parameters (agricultural practices, climatic conditions, soil and vadose zone hydraulic conductivity) can affect surface runoff and infiltration, influencing pesticide transport in soil and unsaturated zone and its fate in surface water and groundwater. Identification of the role of different parameters and evaluation of their influence on the numerical simulation results represent the methodological goal of the study. In addition, the achievement of this goal is important for water management, making possible recommendations for integrated water-pesticide-crop management in the study area.
2. Materials and Methods
Derived from the physically-based analytical model of Hantush and Mariño [38], the adopted mathematical model, with a new Windows-based interface (IPTM-CS), is a hybrid time-continuous and space-discrete semidiscrete model, which takes advantage of both analytical and numerical methods [27,37]. The model is able to deal with physical and biochemical processes related to one-dimensional vertical water flow and three-phase pesticide transport (dissolved, volatilized, adsorbed) in the vadose zone (separated in a root zone and in a deep vadose zone), where complete mixing is assumed. The adoption of a semidiscrete solution allows the simulation for variably-saturated porous media with timeand space-variant parameters under conditions of heterogeneous media, unsteady flow fields, and space-time-dependent physical and biochemical processes concerning pesticide environmental fate. A complete description of the structure of the model and various solution methods are contained in [27,37]. The model has been used in previous studies with consistent results [27,39].
The model is able to simulate pesticide transport both in plant canopy and vadose zone systems, considering different methods of irrigation (overand under-canopy), surface runoff, and possibly of erosion. A one-dimensional, physically-based, compartmental model is used for simulating water flow along the soil profile, divided into surface zone, plant root zone and deep vadose zone, as schematically shown in Figure 2. The new Windowsbased interface of the IPTM-CS model [27] allows easier input and data processing. The large number of tables and figures included provides specific values/ranges of the different input parameters; the reference data are adopted for lacking data or for matching with fieldmeasured data.
In this study, several simulations on representative soil columns were conducted for different situations, starting from real data. First of all, the soil profile was discretized as follows: a) the surface layer was represented by one cell, with an assumed thickness of 0.005 m; b) the root zone was divided in 10 cells of equal thickness, considering a total thickness ranging from 30 to 70 cm for different root zones; as a result, each cell is 1/10 of the assumed total thickness (corresponding to 3 cm for cell when the total thickness is 30 cm, 4 cm for cell in a root zone of 40 cm, etc.); and c) the vadose zone was divided in 10 cells of 20 cm thickness each, for a total of 2 m.
This scheme allowed simulation of different thicknesses of the root zone (from 30 to 70 cm, as field data show), fixing the same structure of the model (10 cells for the root zone), simplifying data input and output. The vadose zone was chosen to have initially a thickness of 2 m, which can be modified taking into account a water table depth ranging from 0.5 to 2.5 m below ground surface. Considering the value which corresponds to the water table depth in the simulated vadose zone, the evaluation of concentration at every possible water table depth can be obtained within the same simulation.
The studied area is 10000 m2, which corresponds to an agricultural unit in the Fucino Plain; the watershed area is about 131 km2, including the entire ancient lake area. The time step adopted was one day for one entire year. Spatial discretization has been assumed appropriate with respect to the time step, considering that the IPTM-CS model does not take into account the process of molecular diffusion and thus evaluation of Peclet number is not necessary. Otherwise, the mechanical dispersion is considered in the simulation.
All performed simulations started on April 1, when agricultural practices including pesticide application at the start of the cropping season occurred. Daily data required for simulation referred to 1989-2004, for the Fucino Ottomila station, located at the center of the Fucino Plain. By taking into account similar trends observed on monthly long-term data (1921-2004), the 1990-93 threeyear sequence, representative of different climate regimen, was used in all simulations. In fact, rainfall in season 1990-91 was higher than the average (+29%), season 1991-92 can represent the average (-2%), and season 1992-93 the low estimate of rainfall (-33%). Temperature daily data were used to calculate daily potential evapotranspiration, using the method proposed by [40].
Simulations were conducted on potato, which is the most common crop in the Plain, representing 23% of the agricultural surface in the last 15 years. Parameter ranges related to potato crops, obtained by databases and tables included in the IPTM-CS program [41–44] are shown in Table 1. Other crops, less represented in the Plain, were not considered in this paper. The crop cycle of potatoes in the Plain helped to determine the values of related parameters. In detail, because potato is cultivated from May to August, both crop coefficient and canopy interception are relevant in this period, while the plant leaf
Figure 2. Scheme of the IPTM-CS model, showing processes taken into account by the model (after Chu and Mariño, 2004).
Table 1. Parameters used for sandy loam potato crop and Linuron pesticide.
area index is active during plant growth (May, June and July).
Characteristics of soils were determined by sampling of root zone soil and further laboratory analyses. The three more representative soils are sandy loam, clay loam, and silty-clay loam by the ASTM classification [45]. The granulometry, the initial water content and the bulk density were laboratory-determined, whereas the saturated hydraulic conductivity was calculated by constant-head, field permeameter tests. From the analysis of soil sections, it was evident that under the root zone there is a more permeable vadose zone, characterized in the sandy loam root zone by the presence of gravels and by a low percentage of silt. Direct surveys did not show fractures or preferential flowpaths in the root zone. This could have resulted in a lower hydraulic conductivity for the root zone (e.g., for sandy loam 6 x 10-6 m/s) and a higher one for the vadose zone (up to 6 x 10-4 m/s). The vertical distribution of the saturated hydraulic conductivity, from the root zone to the lower limit of the vadose zone, matched the average value obtained by in-situ permeameter tests for the entire soil thickness (2 x 10-4 m/s). Different values were also used for other soil parameters, resulting in different characteristics for the root zone and the deep vadose zone in these simulations. Figures and tables included by the IPTM-CS program provided values for parameters which have been not measured and provide methods for the determining these parameters [23,46,47]. The soil parameterization was based on Van Genuchten’s relationships [48]. Soil parameters adopted for simulation in a sandy loam soil are shown in Table 1.
Pesticide parameters for Linuron, Dichloran, and Carbaryl were obtained from public and web databases [49] and tables attached to the IPTM-CS program [14,50,51]. Linuron is a preand post-emergence herbicide used to control annual and perennial broadleaf and grassy weeds on both crop and non-crop sites; it is used frequently for potatoes and carrots. Dichloran is a fungicide widely used for the pre-harvest treatment of several fruit and vegetables crops, including potatoes. Carbaryl is a commonly used wide-spectrum carbamate insecticide of large use, which controls over 100 species of insects on several fruits and vegetables [49].
Table 1 shows values of parameters used for Linuron’s simulation. Use was made of results of previous surveys as targets for the simulated concentrations. In other words, measured pesticide concentrations in 2004 and 2006 in surface water and groundwater were used to determine if the simulated concentrations in soil and in groundwater give similar values. Data of pesticide concentrations obtained by HPLC analysis [33] showed absence of pesticides in detectable amounts in waters during winter (December 2003) and spring (early April 2004). This is in agreement with the agricultural schedule in which the first application of pesticides occurs approximately in mid-April. In summer 2004, Linuron was found in 30 of 35 samples, Dichloran in 6 of 35, and only one case of Carbaryl presence; in Fall 2004, Linuron was detected only in 7 of 35 samples, while Dichloran was largely present in 12 of 35 samples, and Carbaryl was measured in 4 samples [52]. At the end of April 2006, a new survey was conducted in 20 sites selected from the previous 35 sites, where higher concentrations were found in 2004. Linuron was found in 17 of 20 water samples and Dichloran in 7 of 20. Concentrations ranged from 0.03 g/L to 2.8 g/L, with relative abundance in surface waters and lower concentration in groundwater (between 0.03 and 0.5 g/L). June and September 2006 surveys confirmed previous results, having 100% of Linuron occurrence in June (average: 1.72 g/L), while in September less than half off samples have Linuron, with a lower average of 0.1 g/L. Dichloran was found in only one sample in June, having major occurrence in September 2006. In groundwater, mainly Linuron has been found, having the occurrence of 80% in June 2004 and 2006. Occasionally, Dichloran was found in traces (< 0.1 g/L). Samples without trace of Linuron corresponded to groundwater related to deep circulation (deep wells), while in shallow groundwater (shallow wells and springs) Linuron normally was detected, in September 2006 too. Table 2 gives concentration values measured in 2004-06 surveys for pesticides in surface waters and groundwater.
Given that the watershed under consideration is flat, simulations did not take erosion into account. Curve numbers, potential retention, and soil cover factor were evaluated by tables attached to the IPTM-CS program for potato crops and above-mentioned soils [53].
Irrigation practice information was obtained directly by farmers’ interviews and by the agricultural schedule of the Plain. For potato, under-canopy irrigation was performed four times between June and July with an interval of 10 days. The irrigation rate was 0.367 cm/hour, with each irrigation occurring for 4 hours resulting in an application of 1.47 cm.
Finally, pesticide applications, determined by the previously mentioned farmer interviews and by the total amount of pesticides sold in the studied area, were considered as follows: for Linuron, one application at the beginning of April with a rate of 0.117 g/m2; for Dichloran, two applications every 10 days after June 1, with a rate of 0.1 g/m2; and for Carbaryl, three applications every 10 days starting from June 20, with an increasing rate ranging from 0.05 to 0.07 g/cm2. The pesticide application was considered instantaneous and under-canopy because it was spread mainly in non-growing periods. The initial pesticide concentration entered in the simulation files was zero for the first year; for the following years, the value of initial concentration corresponded to the amount obtained by simulation in the last day of the previous year, to consider possible accumulation processes into the soil.
3. Results
Main parameters influencing dissolved pesticide fate and its concentration in waters have been identified as pesticide characteristics, including amount, rate, and single properties of each compound used, recharge (amount and distribution, due to both rainfall and irrigation practice), and physical soil characteristics (e.g. size, saturated hy-
Table 2. Measured pesticide concentrations in water samples (35 samples in 2004; 20 samples in 2006; n.f.= not found).
Table 3. Comparison between simulated and measured Linuron concentrations in surface waters. Simulated values have been obtained considering the average runoff in 15 days (one week before and one week after sampling) (modified from Pacioni et al., 2007). Data in italics highlight periods showing significant discrepancies (April samplings).
draulic conductivity, thickness of root and vadose zones).
To calibrate and validate the model, preliminary simulations of Linuron contribution to the runoff were compared with its pesticide concentrations found in surface waters during all 2004 and 2006 surveys (Table 3). Limited occurrence of pesticide in groundwater did not allow a similar comparison for groundwater. Contaminant dilution in the canal discharge, variability of the applied pesticide amount, and interference of more than one rainfall event have been considered [54]. To take into account possible shifts in the real pesticide application by farmers and the time transit of runoff waters into the Plain, an average of 15 days of the pesticide runoff obtained by simulation has been adopted, considering 7 days before and after sampling. Pesticide amount really used by farmers can not be exactly evaluated in detail for each agricultural unit. These uncertainties certainly affect the comparison between real and simulated concentration, but at the entire basin scale (130 km2) application and runoff processes can be compared, looking for the presence of pesticides in surface waters in a week after the presumed date of pesticide application. Results of this comparison were oriented to highlight the role of runoff, while the specific aim of the research was to simulate pesticide fate in soil and groundwater.
Simulation results showed a significant agreement between simulated pesticide concentration due to runoff and measured concentration found in surface waters. Taking into account the difficulty to verify days and amount of pesticide application, which can locally change with respect to the assumed standard, simulation results were relatively consistent with the monitored ones. Exceptions were noted in April surveys, when the local effect of pesticide application can highly affect pesticide content in surface waters, causing outlier peaks in monitored surface waters (Table 3). Discrepancies between simulated and real runoff data did not affect the reliability of simulations of infiltration to groundwater, as confirmed by specific column experiments [55-56].
The model calibration using runoff data can be accepted for further simulations, taking into account the aim of this research which is to provide a preliminary evaluation of groundwater vulnerability to pollution by pesticide at the basin scale. Considering the adopted methodology as suitable for the study area, several simulations were conducted to determine the effects of hydrological processes on pesticide fate. Firstly, simulations using three different soils (sandy loam, silty-clay loam, and clay loam) and for three detected pesticides were performed, to identify situations that can favor groundwater pollution. In soils with low hydraulic conductivity, like silty-clay and clay loams, pesticides reached the vadose zone only in undetectable concentrations, lower than 0.001 g/L at 0.5 m below ground surface. In contrast, results for pesticide applications in different soils indicated that the runoff process is relevant during irrigation practice and natural decay is largely prevalent for pesticides having low persistence such as Dichloran and Carbaryl (half-life of few days for horticultural crops). For both of those pesticides, simulated concentrations at the base of the root zone were of the same order of magnitude (less than 0.001 g/L at 0.5 m below ground surface) and several orders of magnitude smaller at lower depths. Conversely, application of Linuron, which showed longer persistence in the system, caused highest concentration values in soil pore water. Simulation results (Figure 3) showed the variation of pesticide content with depth in the unsaturated zone. In
Figure 3. Effect of depth on dissolved Linuron concentration with time in a sandy loam soil (year 1990-91). The concentrations shown are those at various depths below ground for a root zone thickness of 30 cm. Simulation period starts on April 1.
the case of high groundwater vulnerability, as for a thin root zone (30 cm) in a sandy loam soil, dissolved Linuron reached 0.1 g/L at 0.5 m below ground and about 0.01 g/L at 1 m below ground (days 260-365 in Figure 3).
Considering this vulnerable scenario, additional simulations were performed using different values and distribution of the abovementioned parameters, starting from a standard scheme that simulates on a potato field the fate of Linuron, the largest occurring pesticide, in a sandy loam soil with 70 cm of root zone and 200 cm of vadose zone.
To evaluate the influence of recharge, different distributions of effective rainfall (rainfall minus evapotranspiration) and irrigation were considered. When it rains in the days following the pesticide application on the ground surface, runoff can easily remove pesticide from the surface and transport it in surface waters. Simulations were conducted assuming non-rainy periods in the following 7, 15, and 30 days since the application of the pesticide (Figure 4). Results indicated that if there is no rainfall after pesticide application, no substantial changes occur in the pesticide concentration in the soil pore water, throughout the year. In fact, the pesticide concentration trend significantly changed only in the spring and summer (days 30-180 in Figure 4), when lower concentration values were registered in case of scarce rainfall. Similar results were obtained considering different years and pesticides. However, this situation which is typical of dry seasons is usually compensated by a major amount of irrigation.
Figure 5 shows a different trend in pesticide concentration obtained from application of higher irrigation amounts. Using as input a double irrigation time, the mobilization of the pesticide into the soil pore water occurred before the autumn infiltration recharge (day 200 in Figure 5), starting in the summer (day 130 in Figure 5).