Pesticide fate in the environment has long been an issue of concern for the civil society around the world. Early ecologist literature, including famous Rachel Carson’s Silent spring (Carson, 1962), has flagged the presence of these chemicals even in areas remote from their application as a strong evidence of the capability of mankind to poison the global environment. Besides the decline in the abundance and diversity of terrestrial organisms, such as Lepidoptera (butterflies), Hymenoptera (honeybees, bumblebees) and beetles (ladybirds, carabids) mainly found in agricultural areas, also freshwater organisms show a high proportion of threatened species as a consequence of exposure to a number of stressors including pesticides (Sala et al., 2000; Vörösmarty et al., 2010). The concern has fostered the ongoing development of stringent legislation on the management of pesticides (EU, 2009b; EU, 2005; European Commission, 2006a; European Commission, 2002; EU, 2009a; European Commission, 2022). The use of pesticides, and related risks, are also considered a criterion to assess the environmental performance of agriculture in the EU (European Commission, 2006b), and is addressed as one of the key aspects of concern in the recent European Farm to Fork (European Commission, 2020b) and Biodiversity strategies (European Commission, 2015; European Commission, 2020a). Moreover, a number of pesticides are included in the list of priority substances targeted by the European Water Framework Directive and related directives (EU, 2000; EU, 2006), which aim at an overall reduction of chemical risk in aquatic ecosystems. Pesticides are widely applied in modern agriculture, typically in sequence during the vegetation period, where pesticide products may contain multiple active substances (Auber et al., 2011; Chiaia-Hernandez et al., 2017). Consequently, in agricultural regions pesticides occur as mixtures in soils (Chiaia-Hernandez et al., 2017), ground waters (Gilliom, 2007; Munira et al., 2018) and surface water bodies (Schäfer et al., 2013; Schreiner et al., 2016). The distribution of pesticides in the environment and their cumulative risks must be quantified when assessing the implementation of EU policies in the field.

Assessments of pesticide risk at European scale have been performed in the past several years based on a variety of approaches (e.g., by Schriever et al. (2007), using an empirical indicator of runoff potential; by Tiktak et al. (2004), using a physical modeling approach; by Delbaere and Nieto Serradilla (2004), using a pesticide fate indicator, among others). Risk assessment usually rests on exposure/toxicity ratios, i.e., environmental concentrations divided by threshold values indicating toxic risks with reference to specified receptors (ecosystems, non-target species, humans) (e.g., Gutsche and Rossberg, 1997; Reus and Leendertse, 2000; FOCUS, 2001; Chen et al., 2002; Reus et al., 2002; Padovani et al., 2004; de Zwart, 2005; HAIR project, undated). Concentrations may derive from direct measurements, but these are affected by limitations (e.g., Wolfram et al., 2021). Therefore estimates based on models, ideally using monitored concentration for comparison and validation, in practice are often the only possibility of quantification, especially when addressing a large area.

When concentrations represent the actual conditions in the environment, we can consider the indicators to represent a “true” risk. Conversely, if they reflect a conventional “reasonable worst case” scenario assuming a given pesticide use, they represent a conditional risk. This is typically the subject of assessments for the authorization of pesticide active substances on the market (FOCUS, 2001).

Currently ca. 300 chemical substances are reported in use in the EU as agricultural pesticides (Galimberti et al., 2020). However, little is known about their actual use patterns and emission rates, hence their expected concentrations in the environment.

These knowledge limitations have so far undermined any attempt at assessing the cumulative environmental risk of pesticides at EU level, although there is expanding evidence that the latter may represent an important pressure on aquatic ecosystems (Oliver et al., 2022). In this contribution, we draw a first EU scale cumulative environmental risk assessment of 148 pesticide active substances (AS) in soil and water across Europe, based on recent estimates of pesticide use derived from available data (Udias et al., 2022). In the following sections, we clarify the method used to compute environmental concentrations of pesticide AS, compare our calculations with available monitoring data, present and discuss the cumulative toxicity resulting from computed concentrations. In the final section, we draw conclusions for pesticide assessment and management in the EU including a set of policy recommendations focusing on identified agricultural regions with high cumulative water toxicity, in support of the implementation of EU Strategies and Targets.

Throughout this paper, we refer to the amounts of pesticides released in the environment as “emissions”. We make use of the estimates of emissions by Udias et al. (2022) for 148 pesticide active substances (AS) in Europe, listed in the Supplementary Material S1, Annex 1, together with their physicochemical properties used for modelling. These represent total amounts of each AS applied in agriculture around the reference year 2015, within EU regions (in kg per year). The regional total application was apportioned to land surfaces within each region on the basis of land cover information, as explained in details in the Supplementary Material S1, Annex 2. This resulted in 148 maps of pesticide AS emission, expressed in kg per hectare and year. Consistent with the resolution of environmental data used for subsequent model calculations, we produced these emission maps with a resolution of 1 km² (see Supplementary Material S1, Annex 2).

The mass of pesticides in soil and the loads to the stream network were predicted using a simple steady state box model (Pistocchi, 2010; Pistocchi, 2013). The model describes the soil as a linear continuous stirred tank reactor; we assume an emission of pesticides that is constant in time, although it is well known that pesticides are applied intermittently and often only during a few days in a year. The implications of this assumption will be discussed later. The model accounts for pesticide removal from the soil due to degradation, runoff and leaching, erosion, and volatilization. From steady state mass balance in soils, we compute the loads of pesticide to the stream network through runoff and leaching as well as erosion. Moreover, we consider that a fixed proportion of emissions reaches the stream network after bypassing the soil, to account for direct losses of pesticides (e.g., through wind drift, dripping from distribution equipment and rinsing of chemical tanks). Direct losses of pesticides are intrinsically difficult to model due to their dependence on local management variables. Here we tentatively set direct losses to 1% of emissions. The implications of this assumption are further discussed below.

The loads of pesticide AS coming from agricultural soils through runoff, leaching and erosion, and the assumed fixed proportion of emissions that bypasses the soil ending up directly in water, are then considered as an input to the stream network. For the stream network, we use a steady state plug-flow model accounting for first-order dissipation. In this way it is possible to estimate the concentration of each pesticide AS at each section of the stream network, taking into account all loads coming from its contributing area subject to dissipation along the respective trajectories. Supplementary Material S1, Annex 2 provides a more detailed description of the calculations used for the mass balance in soils and the stream network. The soil model equations are solved in spatially distributed form using map algebra within a geographic information system (GIS) (Pistocchi, 2014). The stream network model equations are solved in a vector format on a network of sub-basins with a spatial resolution varying from a few to a few tens of squared km (Vogt, 2007), similar to the modelling exercise shown in (Pistocchi et al., 2019). Soil properties and water flow in the stream network are described using European-scale datasets further described in the Supplementary Material S1, Annex 2.

The model is deliberately kept very simple. An alternative could have been to use a complex model incorporating an accurate description of chemical fate and transport processes. When applied to predict environmental concentrations from detailed input data including on pesticide use and observations for calibration and verification, complex models promise to yield the most accurate assessment. Application of these models at European scale has been limited so far to assessing conditional risks (e.g., Tiktak et al., 2002; Holman et al., 2004; Centofanti et al., 2008; Hendley et al., 2009; Urionabarrenetxea et al., 2022) capitalizing on the increasing availability of datasets supporting the definition of model scenarios (e.g. Schneider et al., 2007; Blenkinsop et al., 2008; FOOTPRINT, 2008). Assessing “true” risk at EU scale using complex models, though, is not proportionate to the knowledge available in practice. In the absence of detailed data, models entail assumptions (hence uncertainties) that may jeopardize the advantages of an accurate description of processes. Moreover, when model parameters cannot be calibrated on the basis of observations, models of very different structure and complexity are expected to give results of comparable prediction quality, and use of complex models is only advisable after testing simpler screening-level models (Pistocchi, 2013).

The estimates of emissions by Udias et al. (2022), used as input for our calculation, are an extrapolation of reported emissions in selected countries. For a comparison of model results with monitoring data, observed concentrations of pesticides in water were extracted from the IPCHEM platform (see Galimberti et al., 2020 for details). No suitable data could be found on pesticide concentrations in soils to compare with model results.

Once concentrations are estimated for all 148 individual AS, we can refer to two fundamental models to predict the toxicity of chemical mixtures: the model of concentration addition (CA) for chemicals with a similar Mode of Action (MoA) (Loewe and Muischnek, 1926) and the model of independent action (IA) (Bliss, 1939), also called effect addition or response addition, for chemicals with a dissimilar MoA. The CA model posits that the concentrations of the chemicals in a mixture are exchangeable if scaled by an appropriate effect concentration (Belden and Brain, 2018). On the contrary, IA accounts for the probability that multiple chemicals contribute to an effect. Several studies have analysed the performance of CA and IA when used to predict the joint toxicity of pesticides, indicating overall similar performances, although in different cases one approach has shown a better match with observations compared to the other (Backhaus et al., 2004a; Backhaus et al., 2004b; Belden et al., 2007; Cedergreen et al., 2008; Norgaard and Cedergreen, 2010; Schäfer et al., 2013; Cedergreen, 2014). In the light of this evidence, the literature tends to suggest CA as the first step in ecological risk assessment irrespective of the MoA, as it provides a precautionary but not overprotective estimation for pesticide mixtures (Junghans et al., 2006; Backhaus and Faust, 2012; Schell et al., 2018) while more complex approaches may not be applicable widely due to data limitations (e.g. Kim et al., 2018). Another recent review confirmed that synergisms with high deviations from CA based predictions are rare and the use of CA as default approach is recommended, still keeping in mind some exceptions for specific classes of chemicals (Martin et al., 2021). Therefore, CA has become in practice the standard tool for a first tier pesticide risk assessment in China (Chen et al., 2020), the European Union (Frische et al., 2014), and the USA (Belden and Brain, 2018). In this work, we refer only to the CA model and compute a toxicity indicator as the sum of individual AS concentrations divided by their respective effect concentration. For the latter, we use the median of chronic no observed effect concentration (NOEC) species sensitivity distributions (SSD) provided by Posthuma et al., 2019. For chemicals not covered by these SSD, we refer to the 21-days NOEC for aquatic invertebrates reported in the Pesticide Properties Database (PPDB) developed by the University of Hertfordshire (http://sitem.herts.ac.uk/aeru/ppdb/en/index.htm). The assumed effect concentrations are provided in the Supplementary Material S1, Annex 1, for the 148 AS considered here.

Figure 1 shows the ranges of concentrations of all active substances predicted in the stream network (rivers water), together with the assumed effect concentration. All substances appear to be below the effect concentration. Concentrations in soils show a similar pattern, as shown in the Supplementary Material S1, Annex 2. Liquid phase concentrations in soil tend to be lower than in rivers water because river concentrations reflect also the direct loss of 1% of emissions that we assumed to fall directly on water due to spray drift, dripping and other causes, hence not contributing to soil concentrations. Direct losses may be an important, and in some cases a dominant contribution.

To underline the potential importance of direct losses compared to losses through leaching and soil erosion, we also compute the ratio of combined losses through runoff and leaching to emissions, and the one of losses through erosion to emissions, for all 148 AS (Supplementary Material S1, Annex 2), indicating runoff and leaching losses are mostly below 10%, and rather often below 1% of emissions, while losses with erosion rarely exceed 0.1%.

Figure 2 shows the frequency of occurrence of discrepancies (computed as the ratios between observed and computed concentrations) if within a factor 10, between a factor 10 and a factor 100, or above a factor 100, for each of the AS for which monitoring data were available, across all monitoring sites. We can observe that, using estimated emissions, for a majority of substances (53 out of 86) observed and computed concentrations are within a factor 10 at least in 50% of the cases, while 63 out of 86 AS exceed a factor 100 discrepancy in 20% of the cases or less. However, for 10 substances discrepancies are above a factor 100 in more than 50% of the cases. The scatter plot of the central value of the observed and computed concentration ranges (Figure 3) indicates a tendency of the model to underestimate observations, and generally mirrors the pattern of Figure 2.

Many of the cases where discrepancies are larger can be explained with the quantification limits of the observations. For a large proportion of the monitoring data for some AS, true concentrations are likely lower than reported, since measurements below the limit of quantification (LOQ) have been substituted by the value of the LOQ or LOQ/2 in the underlying datasets. For example, the chemical with the worst performance (Benalaxyl-M) has more than 5,000 observations, but all except 4 in the database indicate a constant value of 0.035 ug/L, whereas the calculated values span a range of much lower values. Ametoctradin, a chemical showing mostly discrepancies beyond a factor 100, has 780 observations, all but one at the value of 0.1 ug/L vis-à-vis calculated values ranging over several orders of magnitude, but well below 0.1 ug/L. Observations for all AS suffer from similar limitations, but when concentrations are higher above or around quantification limits the match with calculated values tends to improve. Scatter plots of observations and modelled concentrations for individual substances across measurement sites are provided in the Supplementary Material S1, Annex 3, for the worst-matching AS.

As a benchmark for the model’s skill at matching observed concentrations, we repeated the model calculations with original reported emissions, in the reference countries (Germany, France, Spain, Italy, Ireland, Denmark, Belgium and the Netherlands) for which these were available. A comparison of the calculation results with observations (Supplementary Material S1, Annex 3) does not show clear improvements over the model based on estimated emissions: than the model with reported emissions: for the latter, 46 out of 82 substances are within a factor 10 in 50% of the cases or more, and 53 out of 82 beyond a factor 100 in 20% of the cases or more, while the scatter plot of central values is slightly worse than with estimated emissions.

Figure 4 shows a map of cumulative toxicity in the stream network, computed under a CA model assumption as the sum of concentrations divided by the respective NOEC for all 148 pesticide AS, with estimated emissions as input. The sum of concentrations divided by NOEC represents the amount of “toxic units” (TU) of the mixture. The Supplementary Material S1, Annex 2 also shows a map of the sum of the corresponding losses from soils due to erosion, leaching and runoff and direct losses (1% of emissions) for all 148 pesticide AS, also in this case divided by the respective NOEC. The sum of NOEC-normalized losses can be read as the losses in mass equivalents of a hypothetical substance with unit toxicity. The two maps, shown in Figure 4; Supplementary Material S1 Annex 2 respectively, visualize the estimated distribution of toxicity in rivers and “toxic load” coming from the catchment, hence the parts of the stream network at highest risk and the parts of the catchments contributing the most to pollution. The two maps are rather consistent with each other and highlight a pattern of higher risk in Southern and eastern countries, but with presence of hot spots in various countries. Figure 5 summarizes the distribution of toxicity in the stream network by country, indicating that a cumulative toxicity between 0.1 TU and 1 TU is widespread across the EU, and a cumulative toxicity of 1 TU is also relatively often exceeded. In aggregated terms, 27.4% of the length of the EU’s stream network is above 0.1 TU and 4.3% above 1 TU due to the 148 pesticide AS considered in the model.

The individual maps of concentration normalized by NOEC for the 148 pesticide AS enable also an analysis of which substances contribute the most to overall toxicity. In Figure 6 we plot the frequency of appearance, across the stream network, of each AS among the top 10 contributors to cumulative toxicity. For instance, a frequency of 0.1 would indicate that an AS is among the top 10 contributors in 10% of the stream network. If we consider the whole EU stream network, including also stretches where the cumulative toxicity is low, we tend to identify the most widespread AS. Conversely, if we focus on those parts of the stream network with toxicity above 0.1 TU or 1 TU, we can identify the substances contributing the most to toxicity. For instance, herbicide Glyphosate appears as the 5th most frequent top contributor with reference to the whole stream network, but only as the 37th when looking at the network with toxicity >10 (and similarly above 0.1 and above 1). However, some substances consistently appear among the top contributors irrespective of the toxicity threshold considered. These include deltamethrin, acrinathrin, cypermethrin, chlorpyrifos, tau-fluvalinate, chlorotalonil, which are mostly insecticides.

Using a simple model, we have presented an analysis of the cumulative risk of toxicity for aquatic organisms due to 148 pesticide AS in the EU. We have shown how more than 27% of the EU’s stream network is affected by a pesticide mixture with concentrations resulting in a cumulative toxicity of 0.1 TU or more, and more than 4% above 1 TU. This indicates that pesticide pollution is widespread, and could be a relevant pressure on aquatic ecosystems in Europe. As a first extensive model-based assessment of pesticides at the EU scale, our quantification is laden with uncertainties further addressed below. The uncertainty affecting model estimates of concentrations of each pesticide active substance is at least one order of magnitude, as shown by the comparison with observations presented above. The uncertainty owes to all parameters used in the calculation. Pesticide properties (degradation half-life, soil-water partition coefficient) are usually characterized in controlled experiments, not necessarily accounting for the variability of environmental conditions. However, in many cases concentrations depend critically on the assumed direct losses of pesticides. An additional uncertainty in the assessment of risks is related to the assumed NOEC, which may also be regarded as uncertain within at least one order of magnitude. The two aspects together cause a potentially very large uncertainty on each of the 148 AS considered in this work. The sum of toxicities could suffer from lower uncertainty, if errors were completely random, as underestimations of toxicity for some AS would be canceled out by overestimations for others. A formal uncertainty analysis of the model, anyway beyond the scope of this work, would add limited benefit for the following overarching reasons.

First of all, we focus on 148 pesticide AS out of a total of about 300 chemical AS in current use. Under a concentration addition model, neglecting substances implies underestimating cumulative toxicity.

Furthermore, our assessment of toxicity is based on concentrations in water calculated by the model. These are quite consistent with reported observations in terms of orders of magnitude and ranges, in the majority of the cases. The model tends to underestimate concentrations in comparison with observations (see Figure 2; Figure 3). The correlation between observed and calculated concentrations, though, is rather weak (Figure 3 and Annex 3). This can be explained by uncertainties and bias related to both observed and calculated concentrations. Pesticide concentrations vary significantly in space and time, and the current sampling schemes with regular grab sampling tend to underestimate true concentrations as relevant exposure periods can be missed (Stehle et al., 2013; Oliver et al., 2022). Indeed, the monitoring data used in this exercise contained frequent records below levels of quantification while calculated concentrations were still significant for cumulative toxicity.

Discrepancies between observations and calculations may owe in part to the simplicity of the model: our steady-state calculation apparently does not account for the highly dynamical nature of certain pesticide applications. A systematic bias of our steady state calculation, in particular, is the approximation of emissions as stationary, while a more realistic model should take into account at least the fact that pesticides are released impulsively, once or a few times in a year. In the Supplementary Material S1, Annex 4, we present and discuss a correction factor to account for the effect of an impulsive release of pesticides when using steady state models initially proposed in (Pistocchi, 2010; Pistocchi, 2013; Pistocchi, 2014). In general, however, the simulation of an impulsive release yields lower annual average concentrations than a steady state calculation, which would further exacerbate our underestimation of observations. Therefore, we should assume higher emissions in order to match the available observations. More complex models, describing processes at a finer spatial and temporal resolution, are not expected to show better predictive capacity, because the key fate and transport parameters (degradation rates and, to some extent, partitioning properties) cannot be calibrated under field conditions (see, e.g., Knäbel et al., 2012). Based on the above considerations, we regard our modelled concentrations as a provisionally acceptable generalization of observed concentrations.

Another aspect to consider is the effect concentration used to compute the TU of the mixture. We rest on existing estimates of the NOEC affecting 50% of the species in aquatic ecosystems, according to the species sensitivity distributions (SSD) modelled by Posthuma et al. (Posthuma et al., 2019), and our assessment is therefore conditional to the validity of this threshold. For example, (EFSA, 2013), recommends applying an assessment factor of 3 for toxicity data based on chronic NOEC based SSDs considering the 5th percentile of the distribution, while here the 50th percentile of the distribution was used and no assessment factor was applied, likely leading to an underestimation of risks. Moreover, for substances with missing SSD-based values (Posthuma et al., 2019), we used the 21-day NOEC for aquatic invertebrates. These NOEC values may not be consistent with the other SSD-based values for those AS which are most toxic towards organisms other than invertebrates.

We flag a relatively high risk from chronic exposure when the mixture toxicity calculated from the NOEC exceeds 1 TU, and a lower but somehow significant risk when it is between 0.1 and 1 TU. However, we could not verify any relationship between our calculated mixture toxicity, hence risk, and impacts on aquatic ecosystems. The NOEC may be about one order of magnitude lower than concentrations at which an acute exposure has observable effects on aquatic organisms [e.g., 67], but the relationship between exceedance of NOEC and the risk of ecological effects is still problematic (Smetanová et al., 2014).

Nutrients and hydromorphological degradation often co-occur with toxicants above ecological risk levels in freshwater ecosystems (Schäfer and Piggott, 2018). Although our knowledge of the combined impacts is limited, stressors can interact in different ways (Schäfer and Piggott, 2018): by exerting a combined pressure on ecosystems (Holmstrup et al., 2010; Stampfli et al., 2013; Link et al., 2017), by influencing the sensitivity of organisms towards other stressors (Holmstrup et al., 2010), and by differentially affecting age stages, populations or species within an ecosystem, triggering disruptions in the population dynamics (Bracewell et al., 2019; Van den Brink et al., 2019). Pesticides frequently co-occur (Matthaei et al., 2010; Moschet et al., 2014; Lemm and Feld, 2017; Szöcs et al., 2017) and interact with other stressors (Holmstrup et al., 2010), (Lange et al., 2011; Magbanua et al., 2013; Piggott et al., 2015; Elbrecht et al., 2016; Jackson et al., 2016; Liess et al., 2016; Magbanua et al., 2016; Chará-Serna and Richardson, 2018; Davis et al., 2018; European Environment Agency, 2018; Bray et al., 2019; Chará-Serna et al., 2019; Juvigny-Khenafou et al., 2020), often in a synergistic way (Holmstrup et al., 2010; Jackson et al., 2016; Liess et al., 2016). Maps of pesticide concentrations as discussed here may support the assessment of combinations of stressors on aquatic ecosystems, which so far have struggled to disentangle the effect of chemical pollutants (Vigiak et al., 2021).

This paper presents what, to our knowledge, is the first EU-scale model-based assessment of the cumulative toxicity of pesticide active substances presently authorized on the market, based on a spatially explicit estimation of pesticide use. We have used a simple, spatially explicit model to quantify the cumulative toxicity of 148 pesticide AS in the stream network of the EU. We find a relatively high frequency of exceedance of safe chronic exposure to pesticides, mainly driven by less than 20 of those substances. While our model is theoretically prone to overestimation, available measurements rather suggest an underestimation of concentrations, hence we do not expect to exaggerate the extent and severity of pesticide environmental contamination. Reference to NOEC and the CA model make the assessment precautionary for the goals we have pursued. In this modelling exercise, we have used an improved estimation of emissions compared to a previous exercise that we regard as an “ancestor” (provided for reference in the Supplementary Material S1, Annex 5). However, this has led only to limited improvements in the overall prediction, consistent with the persisting knowledge gaps on actual pesticide use in Europe. In particular, the data currently available at the EU scale do not allow predicting pesticide AS concentrations more accurately than within one order of magnitude. In the future, an accurate characterization of spatial patterns of pesticide emissions remains key to improve model predictions and the overall assessment of pesticides in the framework of current European pesticide legislation. A recently proposed regulation of statistics on agricultural input and output (European Commission, 2021) requires better harmonization and quality of data from 2025 on, possibly enabling more sophisticated modelling in the future. However, pending an improved representation of emissions, the simple approach presented here may still prove sufficient for a first assessment. For the time being, what we present here is a proof of concept for spatially explicit pesticide risk indicators, which could also be presented in aggregated form (see Supplementary Material S1, Annex 6), in support to the monitoring of policy implementation and EU scale pesticide management, compatible with the available knowledge.