Numerical Modeling of Groundwater Pollution by Chlorpyrifos, Bromacil and Terbuthylazine. Application to the Buñol-Cheste Aquifer (Spain)

Chlorpyrifos, Bromacil and Terbuthylazine are commonly used as insecticides and herbicides to control weeds and prevent non-desirable growth of algae, fungi and bacteria in many agricultural applications. Despite their highly negative effects on human health, environmental modeling of these pesticides in the vadose zone until they reach groundwater is still not being conducted on a regular basis. This work shows results obtained by version 5.08 of the Pesticide Root Zone Model (PRZM5) numerical model to simulate the fate and transport of Chlorpyrifos, Bromacil and Terbuthylazine between 2006 and 2018 inside the Buñol-Cheste aquifer in Spain. The model uses a whole set of parameters to solve a modified version of the mass transport equation considering the combined effect of advection, dispersion and reactive transport processes. The simulation process was designed for a set of twelve scenarios considering four application doses for each pesticide. Results show that the maximum concentration value for every scenario exceeds the current Spanish Maximum Concentration Limit (0.1 μg/L). Numerical simulations were able to reproduce concentration observations over time despite the limited amount of available data.


Introduction
According to the United States Environmental Protection Agency (USEPA), a pesticide is "any substance or mixture of substances whose primary purpose is to prevent, destroy, repel, or control a pest" [1]. Over 500 different pesticide formulations are used in agriculture [2], mainly as herbicides and insecticides, to control weeds and invertebrate pests, and thus improve the crops' quality [3]. The groundwater of the Júcar River Basin (JRB) in Spain is located under an intense agricultural exploitation area in which the use of pesticides is very frequent [4]. The predominant crops are citrus, although there are also irrigated areas dedicated to vegetables, as well as rainfed areas where cereals, olives and vines are grown [5]. In this context, pesticides used in agriculture reach groundwater mainly by dragging and leaching, thus being able to pollute the aquifers, reducing their quality and posing an environmental risk and even one to human health [3]. Several recent studies at the JRB have detected the presence of pesticides in both surface [6][7][8] and groundwater [9][10][11]. Pesticides in groundwater have been found over the limits established by current legislation due to an excessive use of them and poorly optimized application methods. In fact, it is estimated that only 0.1% of the pesticides that are applied really exert their effect on pathogens, since they are applied preventively. This surplus is subjected to different processes: photodegradation, assimilation by crops, runoff into surface waters, soil adsorption, chemical and/or biological degradation and infiltration into groundwater [12]. Therefore, it is essential to study the dynamics of pesticides in groundwater, both in terms of degradation and transport processes. In the scientific literature, several models have been developed to predict the transport and fate of pesticides in the environment. Estimations in the vadose zone can be obtained using PESTAN [13], conservative estimates of pesticide concentrations in groundwater are obtained using SCI-GROW [14], general estimations of pesticide fate and transport are obtained using SUTRA [15] and Hydrus [16], the presence of pesticides in surface waters can be analyzed using TOXSWA [17] and pesticide leaching in groundwater can be simulated using PEARL [18].
These mathematical models take into account an important number of physical, chemical and biological processes, as well as pesticide handling practices in the field. These models try to generalize the knowledge of the behavior of pesticides in the analyzed study area, identifying their most important properties, which can be measured in control stations or in laboratories [19]. Moreover, modeling the environmental fate of pesticides has become an important tool for assessing their potential for water contamination [20]. Pesticide models are increasingly used by EU and US authorities to support decisions regarding the approval of pesticide registration.
The use of numerical models has the advantage of being economical and efficient for the evaluation of pesticides while taking into account the large number of relevant aspects of their use in agricultural practices. The selection of the model is justified by the evaluation or the purpose of the study. In many practical situations where not many data are available, using a simple model with fewer parameters is recommended [21]. Numerical modeling to simulating pesticide behavior is an attractive way to perform environmental assessment of agricultural situations [22][23][24][25]. As the first step of the assessment process, it is necessary to identify which pesticide should be used for a specific climate-crop-soil combination, as well as its application rates, so that both crops and the environment are protected.
This work was performed using version 5.08 of the Pesticide Root Zone Model (PRZM5) model [26], included inside the Pesticide Water Calculator (PWC) [27]. PRZM5 is a flexible software that models the fate of pesticides in the environment using relevant local characteristics of climate, soil, hydrology and crop management. It has been developed to estimate pesticide concentrations in groundwater bodies. Output values are obtained in terms of daily, mean and maximum pesticide concentrations, following the regulatory terms established by the USEPA [28].
PRZM5 must be evaluated in terms of the model internal structure, its scientific basics and its capacity to simulate pesticide fate so it can be applied to research. Therefore, the objective of this study was to use PRZM5 to simulate the fate and transport of three pesticides (Chlorpyrifos, Bromacil and Terbuthylazine) which have been identified inside the Buñol-Cheste aquifer of the Júcar River Basin (JRB). To perform this task, a manual calibration process of the pesticide application was conducted, adjusting the value of the parameter "Amount" of the PRZM5 model, as the exact distribution and applied amount of pesticide mass in the crops were unknown.
To achieve these objectives, it is also necessary to have a detailed and synoptic description of the study area, to have a deep knowledge of the physicochemical characteristics of pesticides and to consider actual values of hydrometeorological, hydrogeological and phenological data. The evaluation of pesticide concentrations is carried out by modeling their behavior on a daily time scale and at a local spatial scale.
Results obtained from this study are important as they provide a more in-depth description of the parameters used in the mathematical modeling of pesticides, while also helping to improve confidence in predictions of pesticide concentrations and facilitating model selection for scientific research and groundwater quality management.

Available Tools for the Numerical Modeling of Pesticide Transport in Soil and Groundwater
Some numerical models are available in scientific literature to approach organic chemicals' fate and transport in the non-saturated zone (NSZ). Table 1 shows a comprehensive list of numerical models for pesticide transport analysis with their main features [29]. These models allow computation of the characteristics of pesticide transport through the unsaturated zone until it reaches the aquifer. Some models even propose control and correction measures once the soil or groundwater is contaminated. The use of these models allows the prediction of the mobility and persistence of pesticides, establishment of the potential risks that they induce in health or the environment. Most models are based on the previous knowledge about the irrigation management strategies of a certain crop and the use of specific pesticides and fertilizers with which water management is optimized.
In this work it was decided to use the PRZM5 model due to the advantage it offers over other similar numerical models. PRZM5 in the modeling process accounts for: (i). Data related to climate, soil, hydrology and crop phenology characteristics at local or regional scale. (ii). Data related to the geometric dimensions and physicochemical characteristics of groundwater bodies. (iii). Data related to the physicochemical parameters of the pesticides to be evaluated (e.g., vapor pressure, solubility in water and molecular weight). (iv). Data related to the contaminant fate and transport characteristics (e.g., photolysis, half-life of the pesticide in the soil, foliar degradation and hydrolysis). The model allows selection of the pesticide application date and the corresponding mass applied to the cultivation fields. (v). PRZM5 output data are shareable in regulatory terms by USEPA as daily, mean and maximum pesticide concentration. (vi). Output data of pesticide concentrations in short and long simulation periods.

The Pesticide Root Zone Model (PRZM5)
In this work the Pesticide Water Calculator (PWC) (USEPA, Washington, DC, USA) [27] was used. PWC is a graphical user interface to interact both with PRZM5 and the Variable Volume Water Model (VVWM). While PRZM5 focuses on the pesticide transport in the unsaturated zone of the aquifer, the VVWM allows modeling of the movement and degradation of pesticides reaching a surface water body. The research described below was performed using PRZM5. The model has been developed to simulate the fate and transport of pesticides in the unsaturated zone of the aquifer on a daily scale and in one vertical dimension [26]. The main output of the model is the pesticide concentration in the first meter of the saturated zone of a user-set constant depth unconfined aquifer, just below the phreatic surface [28]. The model integrates the essential physical, chemical, and biological processes that occur during pesticide filtration with the vertical movement of water through the soil [30][31][32]. The results shown were obtained using PRZM5, which considers the following processes [26]: Chemical Application and Foliar Washoff • Chemical Runoff and Vertical Transport in Soil • Chemical Volatilization The mathematical formulations of the different hydraulic and chemical processes that are taken into account by PRZM5 are described below. Water flow in the unsaturated zone is described by Richards [33] as shown in Equation (1): where: K(θ): unsaturated hydraulic conductivity (cm/s) θ: soil water content h: hydraulic head (m) z: vertical coordinate (m) t: time (s) The mass balance equations that account for the different physical-chemical processes suffered by the pesticide through the unsaturated zone are written accounting for the different phases (dissolved, adsorbed and gas) as shown in Equations (2)-(4) [34]: A∆z where: A : transversal section of the soil column (cm 2 ) ∆z : depth (cm) C w : concentration of contaminant dissolved in water (g/cm 3 ) C s : concentration of contaminant in soil (g/g) C g : concentration of contaminant in gas phase (g/cm 3 ) θ : volumetric water content (volume of pore water/total volume of the sample) (cm 3 /cm 3 ) a : volumetric air content in soil (cm 3 /cm 3 ) ρ s : soil density (g/cm 3 ) t: time (days) J D : mass flux due to dispersion and diffusion in the dissolved phase (g/day) J V : mass flux due to advection in the dissolved phase (g/day)  Figure 1 shows the location of the study area, the aquifer of Buñol-Cheste, located on the eastern side of the JRB inside the Valencia Region (Spain). In summer, the average monthly temperature is 22 • C while in winter it is only 6 • C. Average annual precipitation is 350 mm, varying from 280 mm in the southern part of the aquifer to 550 mm in the northern part. In the driest years, the average rainfall is 150 mm/year, while in the wetter years it can reach up to 750 mm/year. ization of the industry near the metropolitan area of Valencia City and the improvement of the main communication routes.
The aquifer is formed by alternating higher-K materials (mainly gravel and sand) and lower-K materials (silt and clay) spatially distributed as shown in Figure 2. Table 3 summarizes the available geological information inside the Buñol-Cheste aquifer for every geological unit. The Buñol-Cheste aquifer has a total area of 542 km 2 and is characterized by fertile soils and Mediterranean climate that favor agricultural activities. Intense agricultural activity has caused high levels of pesticides to be detected both in ground and surface waters. Furthermore, industrial activity in recent decades has increased due to the decentralization of the industry near the metropolitan area of Valencia City and the improvement of the main communication routes.
The Buñol-Cheste geological region is located in the confluence area formed by the foothills of the Iberian and the Betic Mountain Ranges in Eastern Spain. Therefore, groundwater flows mainly from west to east. According to the information obtained from the current JRB hydrological planning [41], the Buñol-Cheste aquifer can be divided into five zones with similar hydrogeological properties, in which saturated hydraulic conductivity (K) varies as shown in Figure 2.
The aquifer is formed by alternating higher-K materials (mainly gravel and sand) and lower-K materials (silt and clay) spatially distributed as shown in Figure 2. Table 3 summarizes the available geological information inside the Buñol-Cheste aquifer for every geological unit.  Gravels, sands, silts and clays Pleistocene-Holocene High-K Table 3 shows that the Buñol-Cheste aquifer presents an important geological complexity with a great variety of structures and geological features that greatly affect its hydrogeological behavior. The complexity of the geological units is the reason why the hydrogeological functioning is not fully known. The main materials capable of storing and transmitting water in the area are carbonates. In almost all cases the aquifer can be considered to be under unconfined conditions. Only in some observation wells the aquifer shows semi-confined behavior. According to the lithostratigraphic column there are a series of carbonate formations capable of developing aquifers. Table 4 shows the details of the available hydrogeological properties of each aquifer section.    Table 3 shows that the Buñol-Cheste aquifer presents an important geological complexity with a great variety of structures and geological features that greatly affect its hydrogeological behavior. The complexity of the geological units is the reason why the hydrogeological functioning is not fully known. The main materials capable of storing and transmitting water in the area are carbonates. In almost all cases the aquifer can be considered to be under unconfined conditions. Only in some observation wells the aquifer shows semi-confined behavior. According to the lithostratigraphic column there are a series of carbonate formations capable of developing aquifers. Table 4 shows the details of the available hydrogeological properties of each aquifer section.

Available Pesticide Observations from the Monitoring Network
Pesticides have been found in high concentrations in the different groundwater bodies that form the JRB. Following current Spanish regulation (Royal Decree 1514/2009), a groundwater body is considered contaminated when a pesticide or by-product exceeds a Maximum Concentration Level (MCL) equal to 0.1 µg/L or when the total concentration of all these compounds is greater than 1 µg/L [42]. According to this criterion, in the Buñol-Cheste Aquifer, non-compliances due to Bromacil, Chlorpyrifos and Terbuthylazine have been observed. The use of Bromacil and Chlorpyrifos has been restricted or prohibited for years, so it is to be expected that their current concentrations in groundwater will decrease with time. However, pesticides are still detected in concentrations lower than 0.1 mg/L, which still makes it necessary to maintain strict control of these chemicals [5].

Available Pesticide Observations from the Monitoring Network
Pesticides have been found in high concentrations in the different groundwater bodies that form the JRB. Following current Spanish regulation (Royal Decree 1514/2009), a groundwater body is considered contaminated when a pesticide or by-product exceeds a Maximum Concentration Level (MCL) equal to 0.1 μg/L or when the total concentration of all these compounds is greater than 1 μg/L [42]. According to this criterion, in the Buñol-Cheste Aquifer, non-compliances due to Bromacil, Chlorpyrifos and Terbuthylazine have been observed. The use of Bromacil and Chlorpyrifos has been restricted or prohibited for years, so it is to be expected that their current concentrations in groundwater will decrease with time. However, pesticides are still detected in concentrations lower than 0.1 mg/L, which still makes it necessary to maintain strict control of these chemicals [5].   Furthermore, Bromacil concentrations were also found in well 08.140.CA003 (San Álvaro) on the same dates. This spatial-temporal coincidence may be due to the fact that Bromacil was used on one-single common application. Bromacil is a product that has been banned for years and its presence may be due to sporadic use, without discarding possible analytical inaccuracies.
Terbuthylazine In order to explain the causes of non-compliance by Terbuthylazine in the groundwater mass, two key points must be addressed. The first one is the location of the monitoring wells. Although most of the surface of the aquifer is covered by forest lands, the observation wells are located in agricultural areas which must be the source of contamination. Despite Terbuthylazine currently not being used by farmers, it was a basic herbicide in citrus cultivation until recent years and is probably still being used. The second aspect to address refers to the physicochemical characteristics of the pesticides. In this work, the Groundwater Ubiquity Score (GUS) [43] was used as a tool to classify pesticides depending on their risk of leaching into groundwater. Terbuthylazine shows a high GUS index (3.07) which indicates that the pesticide actually reaches the aquifer saturated zone.
Chlorpyrifos leaching potential is medium-high (GUS index = 2.57), so its presence in the aquifer is not surprising. However, the fact that it is detected in a single point repeatedly in time suggests the existence of a very localized use.
The other non-compliance situation detected in this aquifer corresponds to Bromacil, which is a highly soluble herbicide, with high mobility (GUS index = 3.44) and moderate persistence. Bromacil was detected frequently and in concentrations of the order of up to 0.2 µg/L in well 08.140.CA142 (Llano de Cuarte). It is striking that the use of this compound in citrus cultivation has been prohibited since 2002, although the last authorizations for its use were given in 2007. The fact that such high concentrations appear between 2012 and 2014 suggests that, despite its prohibition, it has still been used in some sectors. An alternative explanation may be related to the fact that Bromacil is highly soluble, and it shows high persistence. In such a case, its presence should tend to disappear with time.
Therefore, pesticide observations from the monitoring network justify the need for controlling the presence of pesticides in the aquifer, especially in wells 08.140.CA142 and 08.140.CA002 in order to analyze their evolution with time.

Data Sources
Data used in this research were obtained from different sources: (i). Records of the Valencia Water Regional Authority (JRB); (ii). JRB's Automatic Hydrographic Information System (SAIH); (iii). The Spanish National Meteorological Agency (AEMET); (iv). The Atlas of Solar Radiation in Spain using climate data from Satellite Application Facilities (SAF EUMETSAT); (v). Data and cartography provided by JRB.
The climatic variables studied were maximum temperature (T max , • C), minimum temperature (T min , • C), wind speed (v, m/s), relative humidity (RH, %) and precipitation (PP, mm/d). Table 5 shows the application patterns for each pesticide. The number and dates of the applications were obtained from surveys of farmers [44]. The results of these surveys do not provide information about the exact pesticide doses applied so their distribution is not known with certainty. Therefore, as only ranges of pesticide applications are known, their final value and its distributions were calibrated by the modeler after an iterative process that allowed the establishment the value of the PRZM5 parameter "Amount" for each one of the three pesticides. These amounts of pesticide which are provided annually on the area associated with every observation well are shown in Table 5.  Table 6 shows the information sources from which the weather data were obtained for every pesticide. The location of these hydrometeorological stations is shown in Figure 1. The Evaporation Factor was computed following Meyer's formula in Equation (5) [46]: where: E: daily evaporation (mm/d); C: empirical coefficient taken as 0.36; e a : saturation vapor pressure at the water surface (mm hg); e: air vapor pressure (mm hg); v: wind speed (km/h) measured at 7.64 m above the surface. When dealing with citrus crops, the type of irrigation is under canopy. The maximum amount of supplied water is 72 cm/day, and the Soil Irrigation Depth is equal to the full depth of the Root Zone. The unsaturated zone was vertically discretized in seven horizons for the simulation of Chlorpyrifos in well 08.140.CA002 (La Purísima) and also for the simulation of Bromacil and Terbuthylazine in well 08.140.CA142 (Llano de Cuarte). Tables 7 and 8    Variation in soil temperature affects pesticide degradation. To simulate this effect, it is necessary to calibrate the values of Lower Boundary Condition Temperature and the Albedo. The values corresponding to these parameters after a manual calibration process are 12 • C and 0.2, respectively. In this work a single annual crop cycle is assumed. The emergence, maturity, and harvest dates for the simulation of the three pesticides are 1 April, 1 July and 1 October, respectively. Table 9 summarizes the information included in the model regarding the characteristics of the vegetation cover. For all the three pesticides, it is assumed that the pesticide remains on the leaves after harvest. Afterwards, this residue is mobilized through the unsaturated zone subjected to decay and washing processes.

Results
The use of the PRZM5 model allowed the simulation of pesticide concentrations in the Buñol-Cheste Aquifer accounting for the available data. Once the model was calibrated, the results obtained from these simulations are shown below for each one of the pesticides under analysis. Figure 4 shows the effect of increasing the application of Chlorpyrifos from 0.3 kg/ha to 1.5 kg/ha. It was observed that the evolution over time of the experimental measurements and the simulated values followed the same pattern. As shown in Figure Figure 4 shows the effect of increasing the application of Chlorpyrifos from 0.3 kg/ha to 1.5 kg/ha. It was observed that the evolution over time of the experimental measurements and the simulated values followed the same pattern. As shown in Figure 4, in the period 2006 to 2010 there is an increase in the concentration of Chlorpyrifos, exceeding the MCL considered by Spanish legislation. The results show that, by the end of 2011 Chlorpyrifos breakthrough curves abruptly increased until they reached maximum values. In 2012, concentrations continued to oscillate around 0.10 μg/L while the application continues. Finally, during the years 2013 and 2014 Chlorpyrifos concentrations decrease down to zero. It was clearly appreciated that, in 2014, a decrease in concentration began. Figure 4 also shows that the model has been properly calibrated as the observed values are reproduced by the simulation results. However, not all observations are adjusted in the same way. Some correspond to local maximum concentration points while others are included inside the general trend of the simulated results. The results obtained for the calibrated model show that an annual application equal to 0.94 kg/ha of Chlorpyrifos is the one that provides the best adjustments to the observed data. Figure 5   The results obtained for the calibrated model show that an annual application equal to 0.94 kg/ha of Chlorpyrifos is the one that provides the best adjustments to the observed data. Figure 5 Figure 6 shows the temporal evolution of the Bromacil concentration in Well 08.140.CA142 (Llano de Cuarte) for the period 2008-2018. It can be seen that the concentration curves follow a similar pattern for every type of application, since they all show a maximum value and a symmetric distribution.

Bromacil Simulations
Four different sets of scenarios were designed, considering different annual Bromacil application values (0.50 kg/ha, 0.60 kg/ha, 0.645 kg/ha and 0.70 kg/ha). Bromacil was detected for the first time on 11 April 2009 and disappeared on 6 August 2018, so its presence in groundwater was equal to 2598 days. It was observed that the simulated concentrations exceed the MCL threshold of 0.1 μg/L, showing a significant prevalence throughout all the simulation periods.   Figure 6 shows the temporal evolution of the Bromacil concentration in Well 08.140.CA142 (Llano de Cuarte) for the period 2008-2018. It can be seen that the concentration curves follow a similar pattern for every type of application, since they all show a maximum value and a symmetric distribution.  Figure 6 shows the temporal evolution of the Bromacil concentration in Well 08.140.CA142 (Llano de Cuarte) for the period 2008-2018. It can be seen that the concentration curves follow a similar pattern for every type of application, since they all show a maximum value and a symmetric distribution.

Bromacil Simulations
Four different sets of scenarios were designed, considering different annual Bromacil application values (0.50 kg/ha, 0.60 kg/ha, 0.645 kg/ha and 0.70 kg/ha). Bromacil was detected for the first time on 11 April 2009 and disappeared on 6 August 2018, so its presence in groundwater was equal to 2598 days. It was observed that the simulated concentrations exceed the MCL threshold of 0.1 μg/L, showing a significant prevalence throughout all the simulation periods.  The best calibration of the model was obtained considering an annual Bromacil application equal to 0.645 kg/ha. Figure 7 shows the results obtained for this scenario, for which all the observed values but one are well reproduced. The only observation that is not properly reproduced is the one taken on 3 July 2012. On this date the observed concentration was equal to 0.06 µg/L, while the simulation result overestimates this value (0.20 µg/L). There are many reasons that can cause these differences, including errors in the measurement process. The best calibration of the model was obtained considering an annual Bromacil application equal to 0.645 kg/ha. Figure 7 shows the results obtained for this scenario, for which all the observed values but one are well reproduced. The only observation that is not properly reproduced is the one taken on 3 July 2012. On this date the observed concentration was equal to 0.06 μg/L, while the simulation result overestimates this value (0.20 μg/L). There are many reasons that can cause these differences, including errors in the measurement process.

Terbuthylazine Simulations
Terbuthylazine simulations were performed at Well 08.140.CA142 (Llano de Cuarte) following the same procedure as explained before for Chlorpyrifos and Bromacil. Despite results for Terbuthylazine not being as accurate as those obtained for Chlorpyrifos and Bromacil, the modeling process provides results of great importance to understand the behavior of Terbuthylazine in the Buñol-Cheste Aquifer.

Terbuthylazine Simulations
Terbuthylazine simulations were performed at Well 08.140.CA142 (Llano de Cuarte) following the same procedure as explained before for Chlorpyrifos and Bromacil. Despite results for Terbuthylazine not being as accurate as those obtained for Chlorpyrifos and Bromacil, the modeling process provides results of great importance to understand the behavior of Terbuthylazine in the Buñol-Cheste Aquifer.  The best fit of the observed data is obtained with an application value of 0.20 kg/ha. Figure 9 shows the simulated values of Terbuthylazine concentrations in groundwater for this scenario during the simulation period 2008-2018, after the calibration process had been performed. The maximum simulated concentration value is 0.20 μg/L, obtained on 19 May 2011.
As Terbuthylazine has a high rate of mobilization and persistence, the simulated values oscillate throughout the simulation period in the aquifer.  The best fit of the observed data is obtained with an application value of 0.20 kg/ha.  The best fit of the observed data is obtained with an application value of 0.20 kg/ha. Figure 9 shows the simulated values of Terbuthylazine concentrations in groundwater for this scenario during the simulation period 2008-2018, after the calibration process had been performed. The maximum simulated concentration value is 0.20 μg/L, obtained on 19 May 2011.
As Terbuthylazine has a high rate of mobilization and persistence, the simulated values oscillate throughout the simulation period in the aquifer.

Discussion
PRZM5 was used to simulate Chlorpyrifos, Bromacil and Terbuthylazine concentrations in groundwater in the Buñol-Cheste aquifer. For each pesticide and for every application pattern, the following variables were computed: (i) Peak concentration value (µg/L); (ii) Average concentration value (µg/L); (iii) Number of days in which the simulated concentration is higher than the MCL (0.10 µg/L); (iv) Date on which the pesticide concentration is reduced to zero, and (v) Number of days in which pesticide concentration in groundwater is greater than zero. Table 10 shows the results summary using these five variables. A set of twelve scenarios (three pesticides and four application doses) was designed to analyze pesticide concentration in groundwater. The results show that the maximum concentration value for every scenario exceeds the Spanish MCL (0.1 µg/L). As expected, simulations results justify that the simulated peak concentrations increased when the application dose increased. However, some differences were found when analyzing the behavior of each pesticide individually. Average concentrations are lower than the MCLs for Chlorpyrifos and Terbuthylazine but exceed the MCL for Bromacil.
The number of days for which pesticide concentration in groundwater is higher than the MCL are around 2500 days for Bromacil, despite variations in the application dose, meaning that Bromacil is a more persistent compound. For Chlorpyrifos and Terbuthylazine the number of days for which pesticide concentration in groundwater is higher than MCLs change a lot with the application dose (less than 550 days for Chlorpyrifos and less than 1000 days for Terbuthylazine). Terbuthylazine is more sensitive to application dose changes (e.g., doubling the dose increases the permanence of the contaminant 64 times). The persistence of pesticide in the observation wells is similar for Bromacil and Terbuthylazine (between 3000 and 3400 days) and lower than 2700 days for Chlorpyrifos. To analyze the evolution in time of the concentration of the simulated pesticides the corresponding box plot diagrams with the main statistics were obtained (Figures 10-12). The main effects that were observed on the JRB area are described below. The number of days for which pesticide concentration in groundwater is higher th the MCL are around 2500 days for Bromacil, despite variations in the application do meaning that Bromacil is a more persistent compound. For Chlorpyrifos and T buthylazine the number of days for which pesticide concentration in groundwater higher than MCLs change a lot with the application dose (less than 550 days for Chlorp ifos and less than 1000 days for Terbuthylazine). Terbuthylazine is more sensitive to plication dose changes (e.g., doubling the dose increases the permanence of the conta nant 64 times). The persistence of pesticide in the observation wells is similar for Broma and Terbuthylazine (between 3000 and 3400 days) and lower than 2700 days for Chlorp ifos. To analyze the evolution in time of the concentration of the simulated pesticides corresponding box plot diagrams with the main statistics were obtained (Figures 10-1 The main effects that were observed on the JRB area are described below.   The number of days for which pesticide concentration in groundwater is higher th the MCL are around 2500 days for Bromacil, despite variations in the application do meaning that Bromacil is a more persistent compound. For Chlorpyrifos and T buthylazine the number of days for which pesticide concentration in groundwate higher than MCLs change a lot with the application dose (less than 550 days for Chlorp ifos and less than 1000 days for Terbuthylazine). Terbuthylazine is more sensitive to plication dose changes (e.g., doubling the dose increases the permanence of the conta nant 64 times). The persistence of pesticide in the observation wells is similar for Brom and Terbuthylazine (between 3000 and 3400 days) and lower than 2700 days for Chlorp ifos. To analyze the evolution in time of the concentration of the simulated pesticides corresponding box plot diagrams with the main statistics were obtained (Figures 10-1 The main effects that were observed on the JRB area are described below.   The simulations carried out with PRZM5 in this study identified that the annu plication dose of the pesticide (parameter "Amount" in PRZM5) is the most sensiti rameter when considering the soil and climate conditions of the Júcar River Basin. T fore, the "Amount" parameter must be carefully calibrated when the exact values pesticide applications are not available.

Conclusions
This work provides a first step towards the use of numerical modeling techniqu the analysis of pesticide behavior in the Buñol-Cheste aquifer (Spain). For the first this original research shows results obtained by the PRZM5 model to simulate th and transport of Chlorpyrifos, Bromacil and Terbuthylazine between 2006 and 2018 different observation wells (Well 08.140.CA002-La Purísima and Well 08.140.CA142de Cuarte) where pesticide concentrations exceed the current Spanish MCL (0.10 were selected as sources of concentration data. A set of twelve different scenarios, co ering four application doses for each one of the three pesticides, were analyzed. One of the main purposes of this study is the generation of information and too organize and facilitate the performance of risk assessments of pesticides in the B The simulations carried out with PRZM5 in this study identified that the annual application dose of the pesticide (parameter "Amount" in PRZM5) is the most sensitive parameter when considering the soil and climate conditions of the Júcar River Basin. Therefore, the "Amount" parameter must be carefully calibrated when the exact values of the pesticide applications are not available.

Conclusions
This work provides a first step towards the use of numerical modeling techniques for the analysis of pesticide behavior in the Buñol-Cheste aquifer (Spain). For the first time, this original research shows results obtained by the PRZM5 model to simulate the fate and transport of Chlorpyrifos, Bromacil and Terbuthylazine between 2006 and 2018. Two different observation wells (Well 08.140.CA002-La Purísima and Well 08.140.CA142-Llano de Cuarte) where pesticide concentrations exceed the current Spanish MCL (0.10 µg/L) were selected as sources of concentration data. A set of twelve different scenarios, considering four application doses for each one of the three pesticides, were analyzed.
One of the main purposes of this study is the generation of information and tools that organize and facilitate the performance of risk assessments of pesticides in the Buñol-Cheste aquifer. The original information that was produced during the development of this work includes: (i). a set of databases, according to the parameters required by an unsaturated contaminant transport model; (ii). a sensitivity analysis of the model parameters, including the pesticide annual application dose (kg/ha); (iii). estimation of the annual concentration values of every pesticide under a risk assessment framework for the first time in the Buñol-Cheste aquifer.
The greatest difficulty encountered when carrying out the simulations was related to the calibration of the model parameters (Curve Number, USLE factors, soil adsorption coefficient, among others). Given the lack of sufficient field data, in this study it was decided to perform a manual calibration of the pesticide application doses, described by the parameter "Amount" in PRZM5, following a trial-and-error process, modifying this parameter individually and analyzing the change of the model results. Once the manual calibration process was performed, it can be seen that the model results reproduce the concentration observations.
Chlorpyrifos In conclusion, simulation results indicate that, despite the lack of field information, the PRZM5 model is a valid tool to study the behavior of pesticides in the Buñol-Cheste aquifer as it provides valuable information to stakeholders and environmental authorities.  Data Availability Statement: Data used in this research were obtained from different sources: (i) Records of the Valencia Water Regional Authority (JRB). Restrictions apply to the availability of these data. Data was obtained JRB and are available from the authors with the permission of JRB; (ii) JRB's Automatic Hydrographic Information System (SAIH). https://www.chj.es/eses/medioambiente/SAIH/Paginas/Inicio.aspx; (iii) The Spanish National Meteorological Agency (AEMET). http://www.aemet.es/es/portada; (iv)The Atlas of Solar Radiation in Spain using climate data from Satellite Application Facilities (SAF EUMETSAT). https://www.aemet.es/documentos/es/ serviciosclimaticos/datosclimatologicos/atlas_radiacion_solar/atlas_de_radiacion_24042012.pdf; (v) Data and cartography provided by JRB. Restrictions apply to the availability of these data. Data was obtained JRB and are available from the authors with the permission of JRB.