Can information on past and near-future weather and field conditions predict the safest pesticide application day?

Pesticides applied on drained fields can be transported through drain-connected macropores to surface waters, as affected by weather and field conditions around spraying time. We investigated whether field conditions and past and near-future weather could be predictors of the ecotoxicological risk of pesticide leaching to surface waters and help identify the safest application day. We developed a methodological approach for one climate-soil setting and two herbicides. The agro-hydrological model Daisy and a 3,000-year series of generated weather were used to simulate a total of 369,326 pesticide application days and their resulting hourly drain concentrations, used in the risk calculation. Recurrent neural networks were trained on the simulated data. The trained meta -models selected the safest application day within a 5-day period with an average accuracy of 60 – 80%. The effective risk reductions were only of 3 – 11% for the investigated climate-soil-pesticide settings. Nevertheless, they represented 46 – 86% of the achievable reductions according to Daisy.


Introduction
Diffuse pesticide pollution in surface waters from cropland has raised environmental concerns, as the occurrence of pesticides in surface waters poses a risk for aquatic organisms and threatens freshwater biodiversity (de Souza et al., 2020;European Environment Agency, 2018;Malaj et al., 2014;Vörösmarty et al., 2010). However, due to their significant role in supporting food safety and security (Cooper and Dobson, 2007), pesticides remain widely used, with stable sales in the European Union (EU) over the past 10 years (Eurostat, 2020).
Subsurface drain lines are a significant entry route for pesticides into the surrounding aquatic environment, with an average measured pesticide load in the drains of 1 % of the applied dose (Brown and van Beinum, 2009). In Northwest Europe, an estimated 34 % of the total agricultural land is artificially drained (Blann et al., 2009), and more than 50 % in countries like the Netherlands and Denmark. Therefore, reducing pesticide leaching to drain lines could help mitigate environmental pollution from agriculture.
Macropores, here defined as soil pores larger than about 0.3 mm in equivalent cylindrical diameter (Jarvis, 2007), form critical transport pathways for pesticides in agricultural fields by creating direct connections between the soil surface and deeper soil layers, through which water and pesticides can rapidly flow (Petersen et al., 2004(Petersen et al., , 2003. Especially, earthworm burrows can form direct flow paths between the soil surface and the drains (Nielsen et al., 2015a(Nielsen et al., , 2015b, although pesticides in drains can also originate from the subsoil where they can be stored for several years due to a slower degradation (McKnight et al., 2015;Sandin et al., 2018). Preferential flow is strongly influenced by soil properties such as clay content, structure, and permeability (Beven and Germann, 2013). It can be triggered under specific soil and weather conditions e.g., when intense rainfall events occur on a soil close to saturation and water is ponding at the soil surface above or near the opening of a continuous macropore (Jarvis, 2007;Rosenbom et al., 2015).
Since pesticide transport through macropores depends on the timing of specific conditions, many years of observations are needed to quantify the risk. This makes field monitoring of pesticide concentrations in drain pipes unmanageable, in particular when investigating multiple experimental treatments, crops, and crop rotations. Modelling has therefore been widely used to simulate pesticide transport under various conditions (Diamantopoulos et al., 2017;Dubus et al., 2003;Larsson and Jarvis, 2000). In the EU registration procedure for plant protection products, models are used to predict environmental concentrations in surface waters (EU directive 91/414 andRegulation (EC) No 1107/ 2009). The models included in the guidance report of FOCUS Surface Water Scenarios account for different pathways of pesticide transport, namely spray drift, surface run-off, and leaching to subsurface drain lines (FOCUS, 2001). The latter is simulated with the MACRO model, which includes macropore flow (Jarvis and Larsbo, 2012).
Using the MACRO model to simulate pesticide leaching to drains under 1,000 combinations of weather and application dates, Lewan et al. (2009) found that a reduction in soil water content, relatively to soil saturation, at application time and both short-and medium-term total precipitation after application (5, 5 to 30 and 30 to 90 days) were the most sensitive predictors of the total amount of pesticide loss to drains and maximum daily average drainage water concentration. Likewise, MACRO simulations of pesticide fate conducted by Nolan et al. (2008) for three different compounds applied on nine different soil types, at five distinct application dates, and under six synthetic weather data series, showed that in addition to soil texture and pesticide properties, total winter rainfall, short-term average temperature and short-term cumulated rainfall (≤91 days after application) strongly influenced pesticide leaching to drains.
Similarly to MACRO, the agro-hydrological model Daisy (Abrahamsen and Hansen, 2000;Hansen et al., 2012aHansen et al., , 2012b simulates water and pesticide transport in both soil matrix and macropores (Abrahamsen, 2011). In Daisy, the term biopore is used as the description of macropores mostly refers to those biologically-formed e.g., earthworm burrows and root channels (Holbak et al., 2021). These biopores can be parameterized as drain-connected based on field observations Petersen et al., 2012;Nielsen et al., 2015aNielsen et al., , 2015b. Daisy has been shown to successfully simulate fast breakthrough of pesticides to drains and pesticide concentration in drain water with the description of drain-connected biopores (Holbak et al., 2022a). In a modelling study of the Danish Pesticide Leaching Assessment Programme (PLAP) fields (Rosenbom et al., 2015), transport through drain-connected biopores was also essential to match measured glyphosate leaching to drains (Hansen et al., 2012a(Hansen et al., , 2012b. In contrast to MACRO, transport processes of water and pesticide can be simulated in two dimensions (2D) in Daisy, enabling explicit treatment of horizontal transport to drains and making it possible to simulate different soil properties at and between the drain lines. In general, drainconnected biopores only occur within a 1-m belt above the drain lines Petersen et al, 2013;Nielsen et al., 2015aNielsen et al., , 2015b, which calls for a 2D-model for a more accurate description of biopore flow.
In addition, Daisy includes a dynamic crop module (Gyldengren et al., 2020) where, among other things, growth and phenology are a function of the weather data such as radiation and temperature. This means that both application windows and canopy interception will be a function of weather, which is useful when using Daisy for risk analysis. Most environmental models used in pesticides studies include a user defined leaf area index (LAI) model at best (Jarvis and Larsbo, 2012). In Daisy, crop growth is characterized by a LAI and a development stage. The latter can be converted into a BBCH value (Gyldengren et al., 2020), useful for describing specific pesticide application windows.
Using Daisy, Rasmussen et al. (2015) showed that intense rainfall events of 1 h (with intensities between 13 and 39 mm h − 1 ) combined with initially wet soil conditions resulted in higher glyphosate leaching to drains on a sandy loam soil, as compared to a precipitation regime including a series of smaller rain events over 3 to 9 h prior to an intense rainfall event (of 13 to 39 mm h − 1 ). Similar results have also been demonstrated in experimental studies (Edwards et al., 1993;Shipitalo et al., 1990). Hence, specific field and weather conditions on and shortly after the day of pesticide application may affect the fate of pesticides in soil differently.
The studies of Lewan et al. (2009), Nolan et al. (2008, and Rasmussen et al. (2015) thus suggest that pesticide leaching to subsurface drain lines could potentially be reduced without reducing the applied dose, but by appropriate timing of the pesticide application and by accounting for near-future field and weather conditions. In a companion study, the potential risk reduction of pesticide leaching to drain lines was investigated based on Daisy simulations by changing the day of spraying based on Daisy simulations (Holbak et al., 2022b, companion paper). The risk was defined as a combination of hazard for organisms living in streams (ecotoxicological threshold) connected to drain pipes and exposure i.e., predicted concentration in a hypothetic stream. With 800 setups based on unique combinations of 20-year generated weather series for 3 different climates, 3 pesticides, 2 crops, 2 seasons, and 800 synthetic soils, it was shown that a risk reduction of up to 21 % could be achieved by selecting the best day for spraying within a 5-day application period. Given the limited amount of effective mitigation strategies available to farmers for pesticide leaching to drains (Reichenberger et al., 2007), this risk reduction appeared as a significant potential. The reduction could be achieved using the Daisy predictions for each investigated application day which were based on weather information for 300 days after application. However, at most 10 days of weather forecast is typically available. Thus, such Daisy predictions cannot be made with undefined long-term future weather. As a step towards creating a more practical tool for farmers, we tested in the present paper a methodological approach enabling the use of information from the field and short-term weather forecast for selecting safe pesticide application days for one specific soil and climate setting.
The objective of this study was to develop a meta-model based on Daisy simulations, which could predict the safest pesticide application day in a short upcoming window and thereby minimize the ecotoxicological risk of pesticide leaching to surface waters through drain pipes, based on weather forecast and field monitoring data. We hypothesized that past and near-future weather and field conditions such as precipitation, soil water content, and crop development stage, contain sufficient information to predict the safest application day with regard to the ecotoxicological risk.
We first describe the Daisy model and the calibration steps, define the risk of pesticide leaching to surface waters, and detail the optimization process of application days. Second, we report the performance of the meta-models, trained with both past and near-future information, in reducing the risk of pesticide leaching to drains for surface water organisms. Finally, we analyze the sensitivity of the meta-models to the near-future information. To the best of our knowledge, no previous work has been published on attempting to develop a methodological approach that could potentially be used for creating a decision support system for farmers.

Overview of the methodological approach
This work accounts for a sandy loam soil (representative of roughly half of Denmark) in an Eastern Denmark climate, and for two herbicides applied to winter wheat. The herbicides were the non-mobile and persistent diflufenican (Difl) and the mobile and non-persistent iodosulfuron-methyl-sodium (Iodo), along with its metabolite metsulfuronmethyl (Mets).
In a first step, a 2D soil column was parameterized in Daisy and calibrated to field observations of drain flow. Second, 3,000-year long Daisy-simulations were run with a generated weather series. In the simulations, the herbicides were applied on each day within their respective application window and tracked individually over 300 days after application. After transforming drain concentrations into stream concentrations, the criterion mTU was defined for each application day as the maximum ratio of hourly stream concentration to the regulatory acceptable concentration (RAC) (EFSA, 2013).
This resulted in a total of 369,326 application days with their corresponding mTU and weather and field conditions around the application days i.e., soil water content, temperature, and crop BBCH, as simulated by Daisy for the three scenarios. These simulated data were used to train machine learning meta-models. For any random day in an application window, meta-models were developed to predict the day with the lowest mTU among the upcoming 5 days based on weather and field conditions simulated for the 5 previous and next days. These weather and field conditions consisted of both past -up to 30 days before -and near-future information -up to 5 days after. This thus required a total of 10 days of forecasted information on weather and field conditions: 5 potential application days and 5 additional days of input information after the last application day. These 10 days conformed to the typical period of weather forecast available. Furthermore, as the application timing is hardly flexible for farmers due to the pressure on crops by pests, the optimal day should be given in a short upcoming application period, such as 5 days. A schematic overview of the material and methods utilized for developing prediction meta-models of the risk of pesticide leaching to surface waters is shown in Fig. 1.

The Daisy model
Daisy is an agro-hydrological model that simulates water flow and solute (including pesticides) transport in soil (Hansen et al., 1991;Hansen et al., 2012aHansen et al., , 2012b, crop growth (Gyldengren et al., 2020), and organic matter turnover (Bruun et al., 2003). Water flow may take place in the soil matrix and in biopores. As described in Mollerup et al. (2014), water flow in the matrix is calculated with a 2-D solution to Richards equation and solute transport with the advection-dispersion equation. Heat balance is simulated with the conduction-convection equation to model temperature in the soil (Abrahamsen and Hansen, 2000).
The 2-D domain has a plane-symmetry, with the z coordinate describing the vertical direction and x the horizontal direction perpendicularly to the drain (Holbak et al., 2022a). Vertical transport in biopores is included as a source-sink term in the solution to the matrix equations. Biopores may be parameterized as matrix-ending or drainconnected in the biopore model described in Holbak et al. (2021). The vertical transport in biopores is considered instantaneous, meaning that the time for water to move from the soil surface to the bottom of the biopores is less than one numerical timestep, by default one hour. To simulate artificially drained fields, a drain pipe can be included as a single numerical node in the soil domain (Holbak et al., 2021;Mollerup et al., 2014).
Pesticide degradation follows a first-order differential equation based on a specific half-life time (DT50) and is affected by soil temperature, depth, and water pressure potential (Hansen et al., 2012a(Hansen et al., , 2012b. In this study, the sorption of pesticides to soil organic matter is described as an instantaneous process using a Freundlich adsorption isotherm. Metabolites can be generated as decomposition products based on a specific formation fraction and are described with their physical and chemical properties (e.g., DT50, sorption) like the parent compounds.
Daisy was demonstrated capable of simulating the leaching of Fig. 1. Schematic overview of the methodological approach for developing risk prediction models. RAC stands for regulatory acceptable concentration, D0 stands for the day at which one plans to spray the pesticide.
glyphosate as observed in a Danish field (Hansen et al., 2012a(Hansen et al., , 2012b, as well as leaching of bentazone and imidacloprid to drain lines in a Dutch field (Holbak et al., 2022a).

Daisy 2-D parameterization and manual-calibration
The 2D soil column in Daisy was parameterized and calibrated with the aim to establish a realistic setting where the presence of soil layers with low permeability could trigger preferential flow. The parameterization followed the soil concept from Holbak et al. (2022b, companion paper), which describes a conceptual model of the soil profile. The column was described with 5 horizons down to 200 cm depth including two layers with reduced saturated hydraulic conductivity, Ks, or high dry bulk density ρ b (Table S1 and Table 1), based on measurements from an experimental field site, Rørrendegård, Denmark Petersen et al., 2016Petersen et al., , 2008Petersen et al., , 1997. The lower boundary was described as a 200-cm deep aquitard as defined in Mollerup et al. (2014). The Ks of the aquitard and the pressure head at its bottom -the aquitard's boundary conditions in addition to its thickness -were manually calibrated to match observations of drain flow at Rørrendegård (Petersen et al., 2004). The calibration was performed with weather data from Taastrup weather station (1997-2008) located less than 1 km away from Rørrendegård. The objective of the calibration was to reach a reasonable fit to the observed drain water flow using the soil description concept from Holbak et al. (2022b, companion paper).
The hydraulic properties of the different soil layers were calculated with the pedotransfer function Hypres (Wösten et al., 1998), which estimates the parameters from the Mualem and Van Genuchten's model (Mualem, 1976;van Genuchten, 1980) based on texture, organic matter content, and bulk density ( Table 1, Table S  The simulated domain covered the area from the soil surface (z = 0 cm) to a depth of 200 cm, and from the drainpipe (x = 0 cm) to midway between two drainpipes (x = 8 m, for a 16-m drain spacing). The soil layer properties were homogeneous along the x-axis. The drain depth was set to 110 cm as observed at Rørrendegård (Hansen et al., 2012a(Hansen et al., , 2012b, with a drain-ditch-horizon (backfilled trench above the drains created when burying the drain pipes) from 33 to 120 cm depth horizontally located from 0 to 50 cm above the drainpipe (Table 1) as reported by Nielsen (2010).
Previous measurements of biopores with brilliant blue and smoke at two Danish field sites, including Rørrendegård, were used to parameterize the biopore model Nielsen et al., 2015aNielsen et al., , 2015bPetersen et al. 2012;Petersen et al. 2013;Petersen et al. 2016). Both drain-connected biopores starting in the vicinity of the drainpipe (0-50 cm horizontally above the drain pipe) and matrix-terminating biopores covering the whole field (Table S 2) were included.
As drain flow was only observed on one third of the field due to an underlying low-permeable subsoil, the experimental observations of drain flow from Petersen et al. (2004) were scaled with a factor 3 to compare with Daisy simulations where the whole field was considered to have a low-permeable subsoil (Table S 1).
Daisy has previously been reported to successfully simulate pesticide leaching to drains at Rørrendegård field site based on measurements of ioxynil and pendimethalin in drain water (Hansen et al., 2012a(Hansen et al., , 2012b.

Weather
In order to investigate a wide range of weather conditions, a 3,000year weather series generated based on weather data measured at Copenhagen Airport in the years 1983-2013 was used in the simulations (Rasmussen et al., 2018). The synthetic weather consists of hourly precipitation (mm h − 1 ), and daily values of air temperature ( • C), atmospheric pressure (Pa), wind speed (m s − 1 ), and global radiation (W m − 2 ). The weather series simulates realistic variations from hour to hour, day to day, and year to year, but does not include any long-term climatic change. The hourly precipitation was generated using the stochastic rainfall model RainSim V3  and the other weather variables by the climatic research unit (CRU) weather generator (Kilsby et al., 2007), based on their statistical relationship with rainfall. For a detailed description of the weather series please refer to Rasmussen (2017).

Crop and pesticide management
Simulations were run with fixed and conventional management practices that involved ploughing and harrowing and no crop residues were left after harvest. Winter wheat cultivation was simulated every year. The crop was harvested at either crop maturation or latest on the 20th August. Ploughing was carried out on the 10th September, 2 days before sowing. Seedbed preparation was done on the day of sowing. Only management operations that could affect the water dynamics in the soil were included in the management description. Pesticides were consistently applied at 8 AM on the application days, and at the maximum recommended dosage in Denmark (120 g ha − 1 for Difl and 10 g ha − 1 for Iodo (SEGES and Landbrug og Fødevarer, 2021a, 2021b)) for both autumn and spring applications. Difl can be applied in both spring and autumn on winter wheat, while Iodo can only be applied in spring. This resulted in three combinations of seasons and pesticides namely Difl-autumn, Difl-spring, and Iodo-spring. The application windows of Difl and Iodo were defined yearly based on the crop development stage (BBCH), which is weather dependent and thus varied between years. The application guidelines with regards to BBCH were retrieved from the manufactured product labels of DFF and Hussar OD (Bayer CropScience 2015a; 2015b) for Difl and Iodo, respectively. According to the labels, Difl should be applied from BBCH 0 to 28. We limited this window from sowing time to the 1st of November in autumn because the weather and soil conditions do not generally allow for spraying after that date. The application window of Iodo was set from BBCH 14 to 32. This yielded window lengths of 10 to 53 days for Difl-autumn, 15 to 60 days for Difl-spring, and 11 to 81 days for Iodo-spring.
The standard factors of pesticide degradation related to soil pressure potential, depth, and temperature from FOCUS surface water scenarios (FOCUS, 2001) were used in Daisy. The canopy wash-off coefficient was set to 0.075 (unitless). The wash-off function in Daisy differs from that of MACRO as described in Abrahamsen (2018). The environmental fate parameters of Difl, Iodo, and Iodo's main metabolite metsulfuronmethyl (Mets) are described in Table 2.

Pesticide application days
Water and pesticide transport to drains and weather conditions were logged hourly, and field conditions were logged daily. Only application days fulfilling the recommendations regarding upcoming rainfall events were considered. With an application time fixed at 8 AM, this corresponded to days with no rain until 9 AM for Difl (Bayer CropScience, Table 1 Physical parameters for the Rørrendegård field site, based on data from Hansen et al. (2010), Nielsen (2010) and Petersen et al. (2016). Clay, silt and sand fractions are defined according to the USDA classification. The surface crust is assigned a Ks of half the value than in the A horizon (Table S1).

Risk of pesticide leaching to surface waters
The simulated hourly pesticide concentration in the drains was transformed into a stream concentration using a simple stream model described in Adriaanse et al. (2017) (Eq. (1) and Fig. 2). This model assumes that a theoretical stream with a pesticide-free base flow and a catchment area of 100 ha flows along a one-ha field, where spraying occurs as well as on 20 % of the catchment simultaneously (Adriaanse et al., 2017).
Where C s (g/L) is the simulated hourly stream concentration, m d (g ha -1 h − 1 ) the hourly mass of pesticide in the drains, A tf (ha) the assumed treated area including both the field adjacent to the stream (A af = 1 ha) and part of the catchment (20 ha). Q drains and Q base (L ha -1 h − 1 ) represent the hourly drain flow and stream base flow, respectively, and A c the catchment area. The base flow was set to 0.0864 m 3 ha -1 h − 1 , equivalent to 3.6 L ha -1 h − 1 , based on local Danish data (Henriksen et al., 2003).
For each application day, the hourly stream concentration over the next 300 days was converted into toxic units (TU, unitless), i.e. divided by the regulatory acceptable concentration (RAC) of the most sensitive organism. The RAC corresponds to the Effective Concentration for 50 % effect (EC50) of each tracked compound, divided by a factor of 10 (EFSA, 2013). This is equal to 0.134 and 0.057 µg/L for iodosulfuronmethyl-sodium and its metabolite metsulfuron-methyl, respectively, based on the EC50 of the aquatic plant Lemna gibba measured on growth rate after 7 days of exposure (EFSA, 2016(EFSA, , 2015. For diflufenican, the RAC is equal to 0.045 µg/L, based on the EC50 of the algae Scenedesmus subspicatus measured on growth rate after 72 h of exposure (EFSA, 2008). Because iodosulfuron-methyl-sodium and its metabolite metsulfuron-methyl have the same mode of action, their mTUs were combined according to the addition principle.
The maximum TU or the maximum sum of TUs (for Iodo and Mets) was retained for each application day and referred to in the following as mTU, an indicator of the risk of pesticide leaching to drains for aquatic plants and algae. Preliminary work tested the effect of optimizing diverse criteria and demonstrated that selecting application days by minimizing mTU also minimizes the duration at which the stream concentration remained above RAC. mTU was thus the variable to be minimized in our study when selecting an application day, in order to minimize the ecotoxicological risk of pesticide leaching to surface waters. Note that the results cannot be compared directly with a FOCUS assessment (FOCUS, 2001) for purposes of product authorization.

Selection strategies and risk reduction potential
Several selection rules of the pesticide application day were tested: the "best" rule, where the absolute best day (D best ) as simulated by Daisy is selected, the "random" rule, where a random day D random is selected, and the "first" and "last" rules, where the first (D first ) or last day (D last ) is consistently selected, within a 5-day period. The first and last rules were considered as simple empirical rules given that it is commonly advised to spray as early as possible in autumn and as late as possible in spring when the drains are less likely to be active (Brown and van Beinum, 2009;Reichenberger et al., 2007). The best day, D best , thus refers to the day with the lowest mTU (mTU best ) according to Daisy using the full set of soil, weather, pesticide, and crop parameters.
Within each simulated application window, a random day was selected as a starting point. Among the 5 days following this starting point, application days D random , D first , and D last were picked with the selection rules random, first, and last, respectively. For each rule, it was assessed whether the selected day was equal to D best . The accuracy i.e., the measure of the frequency at which D best was selected when following a specific rule, was computed over all application windows. By repeating this procedure 5,000 times with different starting points, we determined the distribution of the accuracy.
To evaluate the practical use of the selection rules, the performance according to the achieved risk reduction was calculated. When following a specific selection rule in a random 5-day period, the corresponding Table 2 Environmental fate parameters of iodosulfuron-methyl-sodium (Iodo), metsulfuron-methyl (Mets), and diflufenican (Difl) (EFSA, 2016(EFSA, , 2015(EFSA, , 2008  mTU of the selected day i.e., mTU best , mTU random , mTU first or mTU last was retrieved. The 90th percentile of each type of mTU was computed over all application windows, and the distribution was determined for each selection rule after repetition of the procedure 5,000 times. A potential risk reduction (%) was calculated based on the estimated distribution's mean, μ, of the 90th percentile of mTU random and mTU best (Eq. (2)).
The potential risk reduction thus refers to the maximum reduction in the 90th percentile of mTU potentially achievable by consistently picking the absolute best day as simulated by Daisy as compared to picking days at random. However, using Daisy for predicting D best requires expert knowledge on the soil-plant-atmosphere system and detailed information on the past and near-future field and weather conditions at the investigated site. To create simpler and more practical models for risk prediction under the scenarios Difl-autumn, Difl-spring and Iodospring we developed meta-models, in which only weather and field conditions potentially accessible to farmers are required as input parameters.

Characteristics of the meta-models
Machine learning meta-models were trained on sequences of simulated application days. Given any single day in an application window, the meta-models were trained to predict D best (i.e., the day with the lowest mTU, mTU best ) within a 5-day upcoming period. The meta-models take input features describing weather and field conditions of the past and (optionally) upcoming 5 days and then predict the safest application day within the next 5 days based on the mTU given as labels. To account for the dependency of the labels on successive past and near-future features, recurrent neural networks (RNN) sequence models were trained using the application programming interface Keras in Tensor-Flow v2 (Keras | TensorFlow Core 2021, v2.3.0).
The considered RNNs consist of a feature extraction part and a prediction part. The feature-extraction part takes field and weather features of the past 5 days as input and extracts a set of non-linear features using 5 stacked layers of bidirectional gated recurrent units (GRU) (Cho et al., 2014) each extracting 25 features from the output of its preceding layer. The bidirectionality of the GRU layers means that features are processed both starting from the oldest to the newest application day and from the newest to the oldest. Bidirectionality improves RNN's preservation of temporal information and allows the network to learn temporal relations in both the forward and backward direction (Cho et al., 2014).
The extracted features are then passed to the prediction part of the network, which consists of 5 unidirectional (forward processing) GRU layers that output a single vector of length 5. Each element in the output vector corresponds to one of the next 5 application days. The prediction part may optionally also consider field and weather features of the upcoming 5 days along with the inputs received from the feature extraction part. The output vector is normalized using the softmax activation function to allow its elements to be interpreted as probabilities (the softmax function maps each value of the vector to the range [0, 1] and normalizes the sum of all vector elements to 1). The aim of the model is to assign most probability to the best application day D best . In the following, D ML_best denotes the meta-model's predicted best day.
The networks were trained on one-hot encoded labels for the next 5 days i.e., on a target vector of values of 0 s and 1 s, where 1 was assigned to D best and 0 to the other four days. The cross-entropy loss function was used to train the models. This loss function considers the model's output probability for each class and corresponding true label (0 or 1) to determine a loss value to be minimized by adjusting the parameters of the network.
The training was performed over 600 epochs (400 steps per epochs, batch data of size 16 fed during each step) with the commonly used Adam optimizer (Kingma and Ba, 2014) and a learning rate of 5x10 -4 tuned based on the validation performance. A dropout regularization coefficient of 0.3 was used for all GRU layers in the feature extraction part, as doing so was found helpful to prevent overfitting as measured by training-to-validation set performance gap (Srivastava et al., 2014). The performance was measured with categorical accuracy.
We first developed meta-models with past and near-future input features from both the previous and upcoming 5 days. We subsequently investigated their sensitivity towards the use of forecasted features such as weather forecast, by re-training the meta-models with past input features from the previous 5 days only.

Features and labels
Based on the studies of Nolan et al. (2008), Lewan et al. (2009), andKoestel et al. (2012), a series of features describing both the past and near-future weather and field conditions was selected, as potentially explanatory variables of mTU. The features, given in Table 3, included field conditions at spraying time e.g., soil temperature and crop BBCH, and precipitation patterns before and after application e.g., time since last or next rainfall. These features were determined using simulation outputs from Daisy such as the soil water content, and the generated weather series used in the simulations. Therefore, these features can be considered perfectly accurate.
The near-future features measuring the time to the next rainfall or to a certain rainfall intensity were all capped to 5 days. Pearson correlation coefficients were computed between all features (Table S 3, Table S 4  and Table S 5) and a correlation coefficient greater or equal to 0.8 (or Table 3 Features for meta-model training in the three investigated scenarios. A "" or a series in brackets indicates that the feature was used. lower or equal to − 0.8) was used to discard some of them. In particular, soil water content calculated by Daisy was only considered in the uppermost soil layer (0-5 cm) as this was strongly correlated to soil water content in the rest of the topsoil. Each scenario dataset was subdivided into three subsets. A share of 80 % of the dataset was used for the training process: 64 % as a training subset and 16 % as a validation subset. A test subset was created with the last 20 % of the dataset and disregarded until training was over (Table S  6). Before training, all subsets were standardized using the mean and standard deviation of each feature in the training set. Centering and scaling the features in each subset allowed having homogeneous data ranges regardless of the units, with feature means of 0 and standard deviations of 1.

Meta-model accuracy and performance on risk reduction
To measure the performance of the trained meta-models in reducing the risk of leaching to surface waters, the selection rules best, random, first, and last previously introduced were compared to selecting D ML_best , the safest day predicted by the meta-models. This is referred to in the following as the "ML_best" rule. D ML_best thus corresponds to the day with the lowest mTU predicted by a RNN meta-model i.e., mTU ML_best , using a selection of weather and field features only. Within each application window of the test sets, a random day was selected as a starting point of a 5-day period to estimate the distribution of the accuracy and the distribution of the 90th percentile of mTU ML_best (see 2.6 for details).
The effective risk reduction (%) of the meta-models was calculated based on the estimated distribution's mean, μ, of the 90th percentile of mTU random and MTU ML_best (Eq. (3)).
Effective risk reduction(%) =μ 90 (mTU random ) −μ 90 (mTU ML best ) μ 90 (mTU random ) *100 (3) The performance of the meta-models was subsequently calculated as the ratio of the effective reduction (Eq. (3)) to the potential reduction (Eq. (2)). The effective risk reduction thus refers to the reduction achieved by consistently picking the safest day predicted by a RNN metamodel D ML_best . Hence, the performance estimates the fraction of the potential risk reduction the meta-models can accomplish.
Additionally, the performance of the meta-models was assessed for the cumulated leached mass of pesticide to the drains. This was to explore whether the meta-models were biased towards the criteria they were trained on. Instead of mTU, the cumulated mass of pesticide leached to drains over 300 days after spraying (g ha − 1 ) was retrieved for each day selected with the random and ML_best rules (defined as D ML_best and D random ). The potential and effective risk reductions were calculated based on the 90th percentile of cumulated leached mass to drains. Furthermore, to investigate the applicability of the meta-models outside their respective training domain, we calculated the performance of the two meta-models developed for spring application on each other's test sets. This estimated how much the risk of leaching of Iodo applied in spring could be reduced by using the meta-model for Difl-spring, and vice versa. Interchanging meta-models may thus indicate whether a common set of rules could be applied for different pesticides.

Calibration of field site
The Ks of the aquitard and the pressure head at its bottom were manually calibrated to 0.035 cm h − 1 and 280 cm, respectively. With this calibration, the simulated duration of both drain periods matched the observations, while the total flow in the first drain season was underestimated, and in the second season overestimated (Fig. 3). This showed that the soil concept from Holbak et al. (2022b, companion paper) was strong enough to obtain a reasonable fit to drain flow measurements by only calibrating the lower boundary. During the simulated two-year period, 3 % of the drain flow came through drain-connected biopores.

Analysis of 5-day periods
We first present a single 5-day period for Iodo-spring to illustrate the effect of application timing (day 1 to day 5) on the subsequent risk (Table 4). The effect of the selection rules on the risk is then presented for multiple 5-day periods (Fig. 4).
In the single 5-day period, the maximum sum of TUs, mTU, calculated after transforming the simulated drain concentration into a stream concentration, was highest for the last day D5 with a value of 0.067 (Table 4). This was 6 % higher than that of the best day D1 (mTU = 0.063). For this specific 5-day period, the trained meta-model successfully classified D1 as the safest day (highest probability of 0.51). In addition, no leaching occurred in the spring following application (Figure S 3); the drain activity stopped a few days after application and no precipitation fell within 10 days after D5. The first leaching event took place in November with preferential flow through biopores. This event coincided with an increased precipitation frequency and thus an increase in biopore drainage. Thus, for this specific 5-day period D5 appeared to be the most risky as the metabolite Mets had the highest concentration peak of 0.019 µg/L at the start of November.
When considering multiple 5-day periods from the test set of Iodospring, the cumulated distribution of mTU by following the different selection rules, namely random, first, last, best, and ML_best rules, can  be computed (Fig. 3). In this example, 70 % of the application years had a similar mTU regardless of the selection rule. When picking the 90th percentile of mTU, the best and ML_best rules were more distinguishable from the random, first and last rules. Beyond the 97th percentile, the mTU appeared noisy. Thus, the 90th percentile allowed to consider a conservative value of mTU without considering extreme events. However, it seemed to leave only little room for improvement between selecting a random or the best day in the 5-day periods.

Accuracy and transport pathways
The accuracy of the selection rule ML_best i.e., the frequency at which D ML_best was also D best , was always higher than for the rules random, first and last, with a mean accuracy of 60, 69 and 80 % for the three investigated scenarios Difl-autumn, Difl-spring, and Iodo-spring, respectively (Table 5). When selecting random days, the mean accuracy was always 20 %, as expected (1 day out of 5, data not shown). The accuracy of the first and last rules varied between the scenarios. In Diflautumn, the accuracy was higher for the last rule, implying that there was often a downwards trend in the 5-day periods, with D best being the first day 30 % of the time and the last day 39 % of the time. In Diflspring, the first and last rules both had a mean accuracy of about 35 %, while in Iodo-spring, the first day appeared to be the best 59 % of the time against 16 % for the last day. These results implied that 31 % of the time, the best application day was neither the first nor the last in Diflautumn, and 25 % of the time in Iodo-spring.
These trends were contrary to the general assumption that the risk is higher when the drains are likely to be activated i.e., in late autumn and early spring (Brown and van Beinum, 2009;Kobierska et al., 2020;Lewan et al., 2009). In the specific 5-day period for Iodo-spring previously presented (Table 4 and Figure S 3), the best day was the first one, as Iodo had a longer time to dissipate from the canopy while both Iodo and Mets had a longer time to degrade on the soil surface before the first rainfall, compared to the last (and worst) day. The time to first precipitation event and pesticide degradation on soil surface also appeared to have a substantial effect on the risk in Holbak et al. (2022b, companion paper). In the specific example, the drains became inactive a few days after that 5-day application period and pesticide leaching to drains first occurred when the drain season started in November. The simulation of pesticide transport through biopores in winter after a spring application revealed that, if the drains are inactive at spraying time, pesticide may remain in the matrix until the next drain season, as also observed in field experiments (McKnight et al., 2015;Sandin et al., 2018). Field studies with the dye tracer Brilliant Blue have also shown that dye applied at the surface can be transported in a small number of biopores and accumulate in the backfill close to the drain lines . When the next drain season starts, pesticide located in the matrix can end up in biopores and thereby leach to the drain pipe, as water and solutes can move from matrix to biopores (and conversely) in Daisy (Holbak et al., 2021).
For Difl-autumn, the analysis of a specific 5-day application period (Figure S 4) shows that Difl leached through preferential solute transport in biopores about 2 months after application. In the analyzed example the best day was the last day due to a higher LAI at application time, as calculated by Daisy based on the soil and weather conditions (Table S 8). The foliar DT50 (dissipation) of Difl being 20 times smaller than the soil DT50 (degradation, Table 2), spraying on a crop with an (even slightly) greater LAI, decreased the amount of Difl reaching the soil surface and entering the soil matrix before leaching to the drains. Thus, spraying later in that 5-day period ensured that the pesticide would hit the canopy first and reduced the risk of pesticide leaching to surface waters. Such results indicate that, for pesticides with a significant difference in foliar DT50 and soil DT50, spraying as late as possible might decrease the risk of pesticide leaching to surface water as more pesticide will end on the canopy and dissipate. For less persistent pesticides with equivalent foliar and soil DT50s this might not be the case and applying as early as possible in autumn might be the safest, as it generally gives more time for degradation before the drain season starts (Lewan et al., 2009). A similar effect of the LAI at application time could be observed for Diflspring in the same specific 5-day application period as analyzed for Iodo-spring (Figure S 5 and Table S 9).
The significance of canopy interception for pesticide leaching to drains indicates that the continuous variable LAI might have a higher predictive power of the risk than the discrete variable BBCH that was used as a feature in the meta-models. However, the determination of the LAI being more complex, the data might seldom be accessible by farmers. Furthermore, because the growth of weeds is not modelled in Daisy, the present results indicate that pesticide interception and dissipation from the plant canopy may be underestimated.

Risk reduction based on mTU
For the two spring scenarios Difl-spring and Iodo-spring, the distribution of the 90th percentile mTU for the selection rule ML_best was clearly distinct from that of the selection rules random, first and last, and nearly overlapped that of the rule best (Fig. 5). The highest potential and effective risk reductions were obtained for the two spring scenarios (10 and 9 % for Difl and 13 and 11 % for Iodo, respectively, Table 5). The meta-model developed for Difl-spring had the highest performance on risk reduction (85 %), although its accuracy was not the highest. The performance was poor for the first rule in Difl-spring (-42 %) and Iodospring (-3%) despite accuracies of 35 and 59 %, respectively (Table 5). Selecting the wrong day 41 % of the time for Iodo-spring had thus crucial consequences on the performance.
For the scenario Difl-autumn, the distribution of the 90th percentile mTU ML_best greatly overlapped that of the selection rule last (Fig. 5 and Table 5). This led to a performance of 46 % for ML_best and 35 % for the rule last.
With a potential risk reduction up to 13 % for Iodo-spring, the risk as determined by mTU appeared not to be greatly different over the 5 days of the application periods. In their work for developing functional tools for pesticide risk assessment and management, Dubus et al. (2009) found that monthly shifts in the application day were best to reduce the pesticide loss to drains as modelled by MACRO. However, In Holbak et al. (2022b, companion paper), where 800 setups based on unique combinations of generated weather series and synthetic soils were used in Daisy simulations, the potential risk reduction was also highest (up to 62 % reduction) when selecting the best application days within the whole application windows, which were up to 55 day long. Nonetheless, when selecting the best application days within 5-day application periods, the potential risk reduction was still significant (9 to 21 % according to the crop, pesticide and application season setting). The smaller potential risk reductions in the present study were probably due to the specific soil and weather setting used in the modelling as well as to the optimization criterion mTU.
The potential risk reduction was especially low in Difl-autumn (Table 5) so that the 5 days to pick from generally had a similar risk. In their modelling study, Lewan et al. (2009) reported that the reduction potential in the 95th percentile of the maximum daily average drain concentration was greater in spring than autumn (89 versus 47 % reduction), when following specific mitigation strategies. In particular, the 89 % reduction for spring application was achieved by restricting application to days with a cumulated rainfall lower than 10 mm over the following 5 days. In the light of our study, their results suggest that the choice of the risk criterion i.e., mTU or maximum daily average drain concentration, may highly impact the calculated risk reductions.
In addition, the risk reduction for each investigated rule was calculated based on the mean 90th percentile of mTU. This statistic was shown to leave only little room for improvement between the selection rules random and best (Fig. 4).
Nevertheless, while the potential and effective risk reductions were small in absolute terms, it was shown that RNN meta-models can be used Table 5 Mean accuracy, risk reductions, and performance of the specific selection rules. The potential risk reduction is obtained by the best rule, the effective risk reduction by the ML_best rule. The performance is calculated as the ratio of the potential and effective risk reductions. The standard errors were smaller than 0.1% for the accuracy. Best: selection of the actual best day from Daisy simulations, ML_best: selection of the safest day predicted by the meta-model, First: selection of the first day, Last: selection of the last day. to achieve a significant share of the potential reduction as illustrated by the distinct 90th percentile mTU distributions (Fig. 5).

Specificity of the meta-models to the risk criterion and pesticide type
When using the meta-models trained with both past and near-future features to minimize the cumulated mass of pesticide leached to the drains instead of mTU, the scenario Difl-spring showed the highest performance (86 %, Table 6). This is in line with a Pearson correlation coefficient between the cumulated leached mass of pesticide and mTU higher for Difl-spring than for the other scenarios (0.32, p-value < 0.0001), highlighting that mTU and cumulated mass were not independent ( Figure S 2). This also corroborates the fact that, in all three scenarios, the cumulated leached pesticide below 110 cm depth over 300 days after application appeared to be lowest for the best day or at least equal to the other days (insignificant numbers) (Table S 7, Table S 8, and Table S 9). Hence, the best day as defined by the day with the lowest mTU was generally also characterized by a low cumulated leached mass. This showed that minimizing the mTU criterion did not only minimize the pesticide concentration peaks, but also the following pesticide leaching events of lower concentrations and potentially longer duration. These results are in line with those of Holbak et al. (2022b, companion paper), where the best day was generally also the best in terms of cumulated mass leached below 200 cm, based on simulations with 800 unique combinations of soil, climate, application season, crop and pesticide. Likewise, the findings of Lewan et al. (2009) indicate that the cumulated leached mass and the maximum daily average drain concentration of a moderately sorbing and quickly degrading pesticide are correlated.
When applying the meta-model of Difl-spring on the dataset of Iodospring, the performance in terms of the 90th percentile of mTU reduction was considerably lower than when applied on its own test dataset (20 %, Table 7). For the meta-model of Iodo-spring on the dataset of Difl-spring, the performance was only 11 %. The performance of the meta-model for Difl-spring applied on the dataset of Iodo-spring was thus twofold higher compared to the opposite combination. Since the first and last rules for Difl-spring had a similar accuracy, we hypothesized that the meta-model of Difl-spring was less influenced by the trends in the application periods, and therefore elaborated more robust prediction rules than Iodospring. The meta-model of Difl-spring could thus successfully predict the safest day for Iodo-spring applications 20 % of the time. Yet, this indicates that using the Difl-spring model to predict the safest application day for Iodo-spring was just as effective as choosing a random day. This implies that the developed meta-models were not generalizable. For a performance higher than 20 %, different meta-models should be developed for different types of pesticides and seasons. Alternatively, pesticide properties such as DT50 could be included as features in the training process.
The specificity of the meta-models towards pesticide properties is in agreement with the results from Nolan et al. (2008), who simulated pesticide leaching to drains on five different soils with MACRO. They found that the correlation between pesticide leaching to drains and weather conditions differed between pesticides. For a moderately mobile and very persistent pesticide, the 75th percentile of predicted cumulated pesticide loss to the drains was strongly correlated to longterm (>91 days) rainfall and winter rainfall, regardless of the application window (Pearson coefficients > 0.8, p-values < 0.05). For a mobile and slightly persistent pesticide applied in spring, the 75th percentile was both strongly correlated to long-term rainfall and moderately correlated (Pearson coefficients > 0.4, p-values < 0.05) to short-term rainfall (≤91 days).

The need for accurate forecasted data
When trained with past features only from the 5 previous days, the accuracy of the meta-models was up to 40 % lower (the ML_best rule in Difl-spring) compared to the meta-model trained with both past and future features (Table 8). The performance in terms of the 90th percentile of mTU was highest for the ML_best rule in Difl-autumn (21 %). In both Difl-spring and Iodo-spring, the performance was drastically lower without future features and even appeared negative in Iodospring.
When compared to the first and last rules, the accuracies of the metamodels trained without future weather information were slightly higher, but the performances were lower (Table 8). This implies that, in the absence of weather forecast information, the meta-models used the trend in the 5-day periods to learn to predict the best day. Yet, because the best day is not always the first or last day, the meta-models did not consistently pick those and, due to the lack of forecasted data, made substantial mistakes in predicting the safest day. This strongly affected the performance and consequently made the meta-models unable to improve the risk reduction beyond what was obtained using the last rule in autumn and spring for Difl.
With a smaller reduction in the performance of the meta-model in Difl-autumn (from 46 % to 21 %) compared to the spring scenarios (from 85 and 86 to 7 and − 3%, respectively), the lack of future features appeared to affect the autumn predictions to a lesser extent. This could be explained by a smaller daily variation in the weather and in the surface soil water content in autumn, making near-future information of a lower importance compared to spring. In the sensitivity analysis of Lewan et al. (2009) carried out with MACRO, the maximum daily average concentration in drain flow of a moderately sorbing and easily degrading herbicide was more sensitive to short-term rainfall after application (5 days) in spring compared to autumn.
Overall, our results are in line with the fact that the drain activity is highly sensitive to the soil water content and the upcoming rainfall, as drains are active when the groundwater level rises above the drain line level. Likewise, upcoming rainfall controls when pesticides stored on the canopy or the soil surface enter the soil matrix or biopores. However, upcoming rainfall provided by weather forecast is known to be uncertain. The performance of the meta-models with perfect future features thus represents a best-case scenario. By discarding future features, we aimed to estimate the performance under a worst-case scenario. The actual performance with noisy weather forecast would therefore probably lie in between these two bounds. Our results show that, to achieve a performance>20 %, accounting for weather forecast as well as predicted field conditions such as soil temperature is necessary. However, it remains to be studied to what extent the performances would be affected by noisy near-future input features as if obtained from forecasts in a realworld scenario, but also noisy features in general.

Optimization of the choice of spraying day
The meta-models were trained to predict the safest day in an upcoming 5-day period, so that they did not take into account the absolute Table 6 Risk reductions and performance of the meta-models in minimizing the 90th percentile of cumulated mass leached to drains. The potential risk reduction is obtained by the best rule, the effective risk reduction by the ML_best rule. The performance is calculated as the ratio of the potential and effective risk reductions. The standard errors were smaller than 0.1% for the accuracy. Best: selection of the actual best day from Daisy simulations, ML_best: selection of the safest day predicted by the meta-model.  Table 7 Risk reductions and performance of the meta-models in minimizing the 90th percentile of mTU when interchanged for the spring scenarios. The potential risk reduction is obtained by the best rule, the effective risk reduction by the ML_best rule. The performance is calculated as the ratio of the potential and effective risk reductions. The standard errors were smaller than 0.1% for the accuracy. Best: selection of the actual best day from Daisy simulations, ML_best: selection of the safest day predicted by the meta-model.  values of mTU (one-hot encoded labels). By doing so, spraying may have been optimized for irrelevant concentrations (one day was best but all days were ecotoxicologically safe i.e., below the RAC) or optimized even though all days were risky (all days were above RAC, but one day was the least risky). Future works could thus investigate the possibility to predict the absolute values of mTU with RNN meta-models. This would help assessing whether spraying in a given 5-day period actually represents a risk or whether another mitigation strategy should be preferred by the farmer, such as a shift to another 5-day application period, in the product to use, or in the method to control weeds (mechanical).
To better differentiate application days within 5-day periods, more complex models could be used to estimate the ecotoxicological risk in future works. This can be done, for instance, by considering the occurrence of both high and long-lasting peaks. The event duration with regards to ecotoxicological effects might help characterizing the risk. With the considered pesticides Iodo and Difl, a critical event could be described by both a stream concentration above the RAC threshold and a duration of the event of 7 and 3 days, respectively (see section 2.5). The estimated stream concentration could also be coupled to a mechanistic ecotoxicological model such as toxicokinetic/toxicodynamic (TKTD) effect models, where the processes that influence uptake, internal exposure, and lead to adverse effects on specific organisms are described (EFSA Panel on Plant Protection Products and their Residues (PPR) et al., 2018).

Conclusions and outlook
Based on a total of 369,326 application days simulated with Daisy for one field site, two herbicides, and one weather series, RNN meta-models were developed to predict the safest pesticide application day with regards to the risk of pesticide leaching to surface waters. The safest application day was predicted for an upcoming period of 5 days. The meta-models developed for the scenarios Difl-autumn, Difl-spring, and Iodo-spring were shown capable of achieving 46 to 86 % of the potential risk reduction. In addition, the meta-models could achieve 35 to 86 % of the potential risk reduction when considering the total leached mass of pesticide to the drains instead of the maximum normalized hourly concentration. Surface processes such as pesticide interception by the crop canopy, canopy wash-off, and surface degradation appeared as factors of the risk of pesticide leaching to surface waters. Such processes should therefore be further investigated to improve their description in Daisy, and especially with regards to the presence of above-ground plant residues. Furthermore, the results of the present study give an estimation of the meta-model performance either under perfect conditions, where the past and near-future field and weather conditions are fully available and accurate, or under a worst-case scenario, where only past features are available. To estimate an actual performance the robustness of these models with uncertain input features should be investigated. Inserting synthetic noise could both help increasing the robustness of the metamodels to natural noise (Karpukhin et al., 2019;Sietsma and Dow, 1991) and estimate their performance under more realistic conditions. This could be done by using information on the data uncertainty from weather forecast providers. Input features should also be subjected to a sensitivity analysis to determine the most important factors and simplify the usage of the meta-models.
By inter-changing the meta-models for the scenarios Difl-spring and Iodo-spring, the meta-models were shown to be specific to the type of pesticide applied. Hence, future work should attempt to develop more generalized models. For instance, pesticide leaching to drains of other types of pesticide could be simulated under different climates (including climate change scenarios) and on diverse soil types (such as the 800 scenarios from Holbak et al. (2022b, companion paper)). meta-models could subsequently be trained with additional features related to pesticide, weather, and soil properties, as the latter have been shown usable for predicting macropore flow ). Generalized or not, with the overarching goal being to develop a decision support system for farmers, the meta-models should eventually be validated against experimental data to be assessed as reliable.
For the tested soil, weather, and pesticide scenarios and under the assumptions of the Daisy model, the safest pesticide application day appeared predictable by the developed meta-models using accurate information on past and near-future weather and field conditions.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: [Jeanne Vuaille reports financial support was provided by Bayer Cropscience. Jeanne Vuaille is currently employed at the European Environment Agency (EEA). The text in this article solely expresses the author's own views and not those of the EEA.].

Data availability
The link to a data repository has been provided in the manuscript