Probabilistic modelling of inherent field-level pesticide pollution risk in a small drinking water catchment using spatial Bayesian Belief Networks

Pesticides are contaminants of priority concern that continue to present a significant risk to drinking water quality. While pollution mitigation in catchment systems is considered a cost-effective alternative to costly drinking water treatment, the effectiveness of pollution mitigation measures is uncertain and needs to be able to consider local biophysical, agronomic, and social aspects. We developed a probabilistic decision support tool (DST) based on spatial Bayesian Belief Networks (BBN) 10 that simulates inherent pesticide leaching risk to groundand surface water quality to inform field-level pesticide mitigation strategies in a small drinking water catchment (3.1 km) with limited observational data. The DST accounts for the spatial heterogeneity in soil properties, topographic connectivity, and agronomic practices; temporal variability of climatic and hydrological processes as well as uncertainties related to pesticide properties and the effectiveness of management interventions. The rate of pesticide loss via overland flow and leaching to groundwater and the resulting risk of exceeding a 15 regulatory threshold for drinking water was simulated for five active ingredients. Risk factors included climate and hydrology (temperature, rainfall, evapotranspiration, overland and subsurface flow), soil properties (texture, organic matter content, hydrological properties), topography (slope, distance to surface water/depth to groundwater), land cover and agronomic practices, pesticide properties and usage. The effectiveness of mitigation measures such as delayed timing of pesticide application; 10%, 25% and 50% reduction in application rate; field buffers; and presence/absence of soil pan on risk reduction 20 were evaluated. Sensitivity analysis identified the month of application, land use, presence of buffers, field slope and distance as the most important risk factors, alongside several additional influential variables. Pesticide pollution risk from surface water runoff showed clear spatial variability across the study catchment, while groundwater leaching risk was uniformly low, with the exception of prosulfocarb. Combined interventions of 50% reduced pesticide application rate, management of plough pan, delayed application timing and field buffer installation notably reduced the probability of high-risk from overland runoff and 25 groundwater leaching, with individual measures having a smaller impact. The graphical nature of the BBN facilitated interactive model development and evaluation with stakeholders to build model credibility, while the ability to integrate diverse data sources allowed a dynamic field-scale assessment of ‘critical source areas’ of pesticide pollution in time and space in a data scarce catchment, with explicit representation of uncertainties.


Introduction 30
Diffuse pesticide pollution continues to represent a significant risk to surface and drinking water quality worldwide (Villamizar et al., 2020). European Union legislations (Water Framework Directive (WFD) (European Commission, 2000), and the related Drinking Water Directive (DWD) (European Commission, 1998), and Groundwater Directive (European Commission, 2006)) require that concentration of individual pesticides in drinking water must not exceed 0.1 µg L −1 and the total concentration of all pesticides must be below 0.5 µg L −1 . Article 7 of the WFD promotes a 'prevention-led' approach to DWD compliance,35 prioritising pollution prevention at source rather than costly drinking water treatment. Catchment management schemes are therefore now widely adopted by policy makers and water companies to mitigate diffuse pollution by pesticides (and other pollutants) and to improve the raw water quality prior to treatment. However, the effectiveness of such diffuse pollution mitigation measures is uncertain due to the heterogeneous nature of catchment systems, and hence catchment management needs to be targeted to consider local biophysical, agronomic, and social aspects (Okumah et al., 2018). To select and prioritise 40 cost-effective interventions, it is essential to identify and map 'high risk' areas, often referred to as critical source areas (CSAs), i.e. those areas within a catchment that contribute disproportionately large amounts of pollutants to a given water quality problem (Doody et al., 2012;Reaney et al., 2019).
Modelling approaches are commonly used to identify diffuse pesticide pollution risk areas and to help evaluate the 45 effectiveness of mitigation strategies. While process-based distributed models, such as the Soil and Water Assessment Tool (SWAT), have been widely used to simulate transport, fate and risks of pesticides at catchment scale and to evaluate the effectiveness of interventions (Babaei et al., 2019;Villamizar et al., 2020;Wang et al., 2019), their application is computationally costly and often hindered by lack of monitoring data for model calibration and validation. Therefore, various spatial index models have been developed to evaluate the intrinsic vulnerability and risk from pesticide pollution at a range of 50 scales (Kookana et al., 2005;Stenemo et al., 2007;Worrall and Kolpin, 2003). The simplest index-based methods assign scores and weights to a set of spatially distributed indicators (e.g., soil media, recharge rate, and depth to groundwater and contaminant properties), which are then combined into an overall risk index, typically within a GIS environment. An example of such index-based method is the DRASTIC system (Aller et al., 1985), which has been widely used for groundwater vulnerability mapping and for identifying areas most at risk to pollutant leaching. DRASTIC only considers geological and hydrogeological 55 factors but ignores the specific nature of the contaminant(s), and it is therefore classed as an intrinsic vulnerability method.
Several modifications and methods similar to DRASTIC exist in the literature, many of which aim to provide specific vulnerability maps, where the contaminant source and behaviour are also accounted for (Duttagupta et al., 2020;Nobre et al., 2007;Saha and Alam, 2014). Other studies have used simplified 1D solute transport models to develop indices and rankings of potential pesticide leaching (Gustafson, 1989;Jury et al., 1987;Stenemo et al., 2007), while other methods such as the 60 SCIMAP modelling framework (Lane et al., 2009;Reaney et al., 2011) use digital elevation models to derive spatial patterns of relative potential erosion and hydrological connectivity to identify possible critical source areas for diffuse pollution risk.

Study site
Jersey Island (c. 117 km 2 ) (49.2138 °N, 2.1358 °W), the largest in the English Channel group, comprises a plateau with an elevation 60-120 m above sea level (Robins and Smedley, 1998). The climate on Jersey is oceanic with average temperature 100 ranging from 6 °C in winter to 18 °C in summer and mean annual rainfall around 900 mm. A shallow, fractured bedrock aquifer underlies most of the island with a generally shallow depth to the water table increasing to 10-30 m beneath higher ground.
For the most part, groundwater storage and transport is shallow and within the top 25 m of the saturated rock (Robins and Smedley, 1998). Bedrock is Precambrian to Cambro-Ordovician. The west-central part of the island is mostly underlain by the oldest rocks belonging to the Jersey Shale Formation, while a volcanic formation occupies much of the east. Superficial deposit 105 include wind-blown sand, loess and alluvium .  Historically, water resources across the Island of Jersey have been vulnerable to nitrate and pesticide pollution, particularly during the growing season when concentrations in untreated water can exceed regulatory drinking water quality levels. This 115 can adversely affect the raw water quality within impounding reservoirs, requiring the water company to undertake a series of mitigation measures at the treatment works to avoid breaches in the treated water supply. The small size of the Island means that land-use across the island is dominated by intensive agriculture, primarily potato and dairy farming. The cultivation of the Jersey Royal Potato takes place during the growing season from January to May. There is very little crop rotation and accordingly the potato crop is grown with the support of man-made fertilisers and pesticides (herbicides, fungicides and 120 nematicides).
This study focused on the Val de la Mare (VDLM) catchment (3.1 km 2 ) in south-west of Jersey, which feeds into the VDLM Reservoir, the second largest reservoir in Jersey constructed in 1962. The reservoir holds up to 938,700 m 3 of untreated water, enough to supply Jersey with water for approximately five weeks. Water feeds into the reservoir from within the catchment 125 area, as well as from neighbouring catchments and a desalination plant when it is in operation. The water quality in the VDLM reservoir is vulnerable to pesticide pollution, with pesticide concentrations often exceeding the regulatory drinking water quality levels of 0.1 µg L 1 for individual pesticides and 0.5 µg L -1 for total pesticide concentration (European Commission, 1998).

130
To evaluate the spatial and temporal variability of the intrinsic pesticide pollution risk from critical source areas within the VDLM study catchment, we have developed a probabilistic model based on spatial Bayesian Belief Networks (BBN). The data and information used to inform the model development are presented in the following sections (2.2 -2.4), while the BBN model itself is described in section 2.5. For the purpose of this paper, we define risk as the probability that the pesticide exposure (i.e., the pesticide loads from the fields to the reservoir) results in the regulatory drinking water standards to be 135 exceeded.

Pesticide detection, usage and properties
Jersey Water, the sole water company and provider on Isle of Jersey, regularly tests and scans the quality of the raw water from the reservoir offtake for a wide range of different pesticides. Five active pesticide ingredients currently or recently in use in the catchment showed evidence of significant concentrations in the reservoir offtake for the drinking water supply. These 140 included the herbicides glyphosate, metobromuron, pendimethalin and prosulfocarb, and the nematicide and insecticide ethoprophos. During 2016-2019 the sampling frequency of these pesticides was not the same (Table 1) with pendimethalin and ethoprophos being sampled on approximately a weekly basis all year round, whereas prosulfocarb generally was sampled on a weekly basis from mid-February to end of May, bi-weekly in January and otherwise every 3 to 5 weeks. The monitoring of glyphosate and metobromuron was more sporadic and only done during spring and early summer with no sampling taking 145 place in the autumn and the winter. Metobromuron was most frequently observed above the drinking water standard, followed by ethoprophos, prosulfocarb and pendimethalin (Table 1). Ethoprophos was not included in the final model, as its use has now been discontinued. In consultation with Jersey Water, it was instead decided to include the nematicide fluopyram as it was considered to be a likely replacement for ethoprophos and hence a potential future risk to the reservoir water quality.
Monitoring data for fluopyram are not available as it is currently not widely used in the catchment. However, fluopyram is 150 known to be used at lower application rates than ethoprophos, making its potential risk of contaminating the water supply intrinsically lower, notwithstanding its relatively high mobility and greater persistence.
Data on pesticide application rates and timing for 2016-2018 was obtained from the Jersey Royals Company, the main potato crop grower in the area who manage ca. 50% of the catchment area, alongside agronomic data on crop coverage and crop 155 rotation for the 2010-2018 period. Pesticide usage (Table 2) was estimated assuming that the available agrochemical data was representative of the whole catchment. Glyphosate was mostly applied in January (mean day of application 23 rd January), while other pesticides were typically applied in February (mean day of application 11 th to 16 th February). Hence, only months January -March were represented in the model, to allow for a potential one-month delay in pesticide application. Usage of pesticides on other crops in the catchment was limited to the use of glyphosate for spraying off grass prior to cultivation and 160 use of pendimethalin on barley. Key pesticide properties for assessing the risk to water quality were extracted from a publicly available database (Lewis et al., 2016) (Table 2). These included the Koc coefficient that represents the adsorption of the pesticide onto the organic carbon of soil, subsoil and vadose zone materials and the half-life (DT50) that represents the degradation rates during transport through 170 each of these layers (Table 2). A third process, that of volatilisation was not considered as this is relatively minor in most cases and omission of this process will provide a conservative estimate for impact to groundwater. Retention of pesticides by the soil and subsoil materials to which it is applied depends on adsorption to organic and mineral surfaces. Soil sorption processes are complex and have been the subject of substantial research in the past. Pesticide sorption is influenced by both the chemical characteristics of the pesticide and soil specific properties, such as soil organic carbon (SOC) concentration, clay content, pH 175 and soil moisture content. For neutral non-polar pesticides, it is well documented that pesticide retention is strongly correlated to SOC concentration, with other factors such as clay content, pH and aeration status playing a subsidiary role (Kah and Brown, 2006;Wauchope et al., 2002). For weakly ionisable pesticides with ionic equilibrium constants near the range of soil pH, sorption may be highly sensitive to the pH of the sorbing soil. However, none of the pesticides considered for the modelling here are considered weakly ionisable, and therefore only the role of SOC content on pesticide retention was considered in the 180 model by including the pesticide adsorption coefficient Koc. The degradation rate of pesticides in soil and subsoil depends principally on the microbiologically-mediated decomposition and as such is also strongly influenced by SOC (Jury et al., 1987;Kah and Brown, 2006). Some chemical degradation of pesticides on mineral surfaces can also occur, which may be more important for the subsoil and vadose zone, but in most cases, biological degradation is seen as the main pathway (Fomsgaard, 1995). In most cases, the time for 50% disappearance (DT50) or 90% disappearance (DT90) is determined. 185 Table 2 Summary of pesticide properties considered in the model (Lewis et al. 2016) and mean application rates in the study catchment. AI=Active Ingredient. Koc = pesticide adsorption coefficient on soil organic carbon. DT50 = time for 50% of pesticide to be degraded in soil.

190
Type of pesticide

Catchment characteristics
Monthly rainfall and temperature data were available from a single meteorological station in the study catchment, operated by Jersey Water for the period 2014-2019. Due to the short duration of this record, additional monthly total rainfall and mean monthly temperature data (from 1894-2019) were obtained from the Government of Jersey website (https://opendata.gov.je/organization/weather). Data for the years 1981-2019 were then used in combination with the 195 catchment-based meteorological data to calculate monthly mean and standard deviation as priors for the model. Mean monthly potential evapotranspiration (PET) and actual evapotranspiration (AET) were calculated using the approach in (Pistocchi et al., 2006), see Appendix A.
Spatial environmental data were processed and collated in a single GIS shape file. Visualisation, geographical analysis, and 200 processing of spatial datasets were done in QGIS 3.12. Environmental Systems Research Institute) was used for the generation of the digital terrain model. The following data sets were used to inform model parameterisation.
Land parcels that fell within the VDLM catchment area were selected and filtered using Feature types (i.e. cultivation) to identify cultivated fields. The field selection was supplemented by visual inspection using satellite imagery (Google Satellite service) to ensure that only cultivated fields were selected. This resulted in the selection of 200 fields, which were assigned a dominant crop type for the 2010-2018 period by spatially joining field polygons with layers containing crop type information.
Crop operation information available for 56 fields within the catchment was used to determine the timing and pesticide 210 application rates (Table 2) and to inform the prior distributions in the model.
A 1m resolution hydrologically-corrected DTM of the VDLM catchment area was created using digital line contours at 1 meter interval provided by Jersey Water and the 'Topo to Raster' tool in ArcGIS; the DTM was used to calculate mean elevation and slope (in degrees) for each field polygon within the VDLM catchment. Topographic connectivity was derived by calculating 215 the horizontal distance from the polygon vertex nearest to the stream features using the Distance to nearest hub tool in QGIS.
Because the study catchment is flat with relatively uniform slopes, we opted here to base the topographic connectivity on the simple distance-to-stream approach. This approach has previously been shown to be accurate for assessing connectivity (e.g., Grauso et al., 2018). For other catchments, it might be more appropriate to use a connectivity based on actual flow direction.

220
The overall depth of the soil column is a key characteristic of the soil that influences the attenuation of surface applied pesticides. Depth to groundwater was calculated for each field using a 5 m groundwater level contour map of Jersey prepared by the British Geological Survey and provided by Jersey Water (T. de Feu, pers. comm, May 20, 2020) by importing and georeferencing the map in QGIS, digitising the contour lines and combining them with field polygons.

225
Soil water retention, conductivity, natural drainage, and anthropogenic characteristics such as plough pans were also considered in model parameterisation. There has been no systematic survey of the soils of Jersey in the traditional sense of classifying and grouping soils according to their pedology. Brief descriptions of 'soil series' were given in Jones et al. (1990) and the Soil Atlas of Europe (European Soil Bureau Network. Eds.: Jones et al. 2005) shows a single soil type Dystric Cambisols for the whole Island of Jersey at a 1:2 500 000 scale. Due to this lack of detailed soil mapping, the hydrogeological 230 map of Jersey  was used to identify soil hydrological units based on the three hydrogeological formations identified in the VDLM catchment (Fig. 1b). These three soil hydrological units consisted of soils developed on loess, aeolian (blown) sand and shales, respectively. Both loess and sand are periglacial Quaternary deposits that are relatively common throughout northern Europe. Soil hydrological data for similar soils is contained within the HYPRES soil hydrological database (Wosten, et al. 1999;Wosten et al. 1998), hence this database was used to derive the soil hydrological properties necessary for 235 the modelling of pesticide attenuation.
According to the descriptions of soil series in Jones et al. (1990), the soils developed on these three hydrogeological units are generally well to moderately well drained with no real inhibition to the downward movement of water. Soil hydrological properties derived from HYPRES database supported this assumption, with mean saturated hydraulic conductivity largely 240 greater than 10 cm day -1 . Less than this would indicate the presence of a slowly permeable horizon (MAFF, 1988) with some degree of ponding within the soil.
Local knowledge suggested that some of the soils in the catchment could have a plough pan potentially as deep as 60 cm below the surface. This layer is likely to restrict the downward flow of water through the soil and increase the likelihood of near 245 surface runoff and was included as one of the soil parameters that could be manipulated within the model. In the absence of data in the HYPRES database on plough pans, a value of 0.02 cm day -1 was selected as a worst case estimate of plough pan hydraulic conductivity based on Koszinski et al. (1995) who reported hydraulic conductivity of compacted soil of between 0.02 and 3 cm day -1 .

250
Laboratory measurements of topsoil soil organic matter (SOM) content and pH were available for 40 and 37 fields within the VDLM catchment, respectively, most of which were cultivated with Jersey Royal potatoes (32). Soil organic matter (SOM) content ranged from 1.9% to 3.8% with mean SOM at 2.5%. Mean SOM for loess was 2.5% (1.9-3.8%) and 2.3% for shale (1.9-3.4%), Median SOM was slightly greater for the six fields found on grassland (3%) than in the potato fields (2.3%). Soil pH ranged from 5.3 to 7.2 with a mean pH value of 6.4 and mean soil pH was slightly greater in the potato fields (6.4) than in 255 grassland (6.1). However, as this information was not available for all fields and only for the topsoil, these data were not used in the model parameterisation. Instead, model sensitivity to using observed SOM concentrations vs. those derived from the HYPRES database was examined and was found to be negligible, therefore the converted HYPRES SOC values (by dividing SOM by 1.724) were used as priors in the model to ensure spatial consistency.

Field attributes
The spatial data described in 2.3 was used to inform the parameterisation of the model variables. It was found that most study fields within the VDLM catchment lay in relatively flat ground with 152 fields having a mean slope less than 3 degrees, while elevation range was also relatively small, from 67m to 99m. Loess was the dominant soil parent material in 155 fields, while soils derived from Jersey Shale or Blown Sand were dominant in the remaining 31 and 14 fields, respectively. Loess covered 265 the central and northern part of the VDLM catchment, while Shale was found in the south-western part and Blown Sand in the south-eastern part of the catchment. Most study fields had a sandy silt loam (120) or sandy loam (59) soil texture class. Fields with a sandy silt loam texture were mostly underlain by the Loess hydrogeological formation (101), while fields with sandy loam texture were mostly underlain by the Jersey Shale formation (36) and blown sand (12).

270
The horizontal distance of VDLM fields to the stream network was used to represent field connectivity to the stream network and the reservoir. The median horizontal distance was 119 m, with 91 fields located within 100 m of the stream network. Depth to groundwater within the VDLM catchment ranged from 0 to 20 m, but most of the catchment area (56%) had a shallow aquifer less than 5 m deep (117 fields), with further 54 fields having groundwater depths between 5 and 10m.

Spatial Bayesian belief network risk model
We have developed a probabilistic model, based on spatial Bayesian Belief Networks (BBN) (Appendix A), to evaluate the spatial and temporal variability of the intrinsic pesticide pollution risk from critical source areas within the VDLM study catchment. The model aims to provide a field level assessment of the relative water pollution risk characteristics of each field, made available as probabilistic map layers. The developed approach integrates the various information sources described above 280 and includes causal relationships between both discrete and continuous variables in a hybrid BBN. The general principles and theory of Bayesian networks have been described extensively elsewhere (Korb and Nicholson, 2010;Moe et al., 2021) and will not be discussed in detail here.
The model structure and development were informed by expert knowledge and stakeholder feedback. Hydrologically, the 285 VDLM reservoir was assumed to be fed by the west and east streams ( Fig. 1) as well as by groundwater, and the groundwater aquifer was assumed to be unconfined and homogenous. Thus, pesticides applied to a given field could either leach to groundwater or could be transported directly to the reservoir or one of the streams through surface runoff. The risk assessment model therefore accounted for pesticide losses via both groundwater leaching and overland flow, with the final assessment based on the combination of both. Both transport pathways are influenced by three key factors, namely soil and site 290 characteristics, climate and hydrology, and land management (e.g., pesticide usage and properties, land use) (Fig. 3, Appendix A).
In the following, the modelling of the groundwater leaching and the overland flow risk components is described. In both cases, we are interested in determining how much of the applied pesticide mass L0 [kg ha -1 ] may eventually reach either groundwater 295 via leaching (Lgw [kg ha -1 ]) or surface water via overland flow (Lof [kg ha -1 ]) from a given field. To ensure mass balance, we assumed that for each field only a proportion fleach of the applied pesticide mass would be available to leaching, while the remaining proportion would run off to surface water. During transport to groundwater and surface water the pesticide can undergo attenuation with the degree of attenuation determined by attenuation factors, respectively AFgw (Eq. 10) and AFof (Eq. 11) (as described in the following sections). Overall, the pesticide loads that reaches groundwater via leaching and surface 300 water via runoff from a given field was therefore given by: Conceptually, this way of calculating pollution risk to groundwater and surface water is similar to the pesticide impact rating index proposed by Kookana et al. (2005) and to the InVEST nutrient delivery model (Sharp et al., 2020). For the modelling here, it was assumed that the fraction of the applied pesticide to land that would be available to leaching fleach equalled the ratio of infiltration to excess rainfall. 310 The combined pesticide load from a given field was the sum of the leaching and the overland flow component.
This combined pesticide load was converted to a surface water concentration Csw [µg L -1 ] to evaluate the risk to surface water as follows: where Ac is the total area of all fields in the catchment (192 ha) and Vres is the water volume in the reservoir (938,700 m 3 ).
Equation 3 is a very simplified way of evaluating the field-level risk; it is essentially assumed that the combined load from the field in question represents the average load from all fields in catchment. This allows the combined pesticide load from the 320 given field to be converted to a concentration in the reservoir, which can then be compared to the regulatory standards based on which the risk can subsequently be assessed. Hence, if the combined pesticide load from a field resulted in Csw exceeding the standard of 0.1 µg L -1 , this field was considered high risk (see Appendix A).

Groundwater leaching risk assessment 325
The conceptual model for the pesticide leaching to groundwater applied in this study catchment ( Fig. 3) followed on from the screening model proposed by Jury et al. (1987). The model assumes that a single pesticide mass flux is applied to the soil surface (z=0). The pesticide is assumed to dissolve in the infiltrating rainwater and move downward through the soil profile by leaching at a constant infiltration rate Jw [m d -1 ], which is here determined by the amount of excess rainfall and the physical properties of the soil (Appendix A). The pesticide is assumed to move downward through the soil by piston flow (i.e. no 330 dispersion) while undergoing linear adsorption and first-order decay. Given these assumptions, the pesticide transport and fate can be described by the following mass balance equation (Jury et al. 1987): where C [µg L -1 ] is the pesticide concentration in solution (the water phase), w is the volumetric water of the soil, µ [d -1 ] is the first-order degradation constant, and RF is the retardation factor given by Jury et al. (1987): where Koc [L kg -1 ] is the organic carbon partition coefficient, and b [kg L -1 ] and foc [kg kg -1 ] are the soil bulk density and the organic carbon content of the soil, respectively. The retardation factor describes the velocity of the solute pesticide relative to the infiltrating water. The solution to the mass balance equation (Eq. 4) for an instantaneous mass injection m0 [kg ha -1 ] can be written as: where m(z) [kg ha -1 ] is the pesticide mass contained within the pulse when it reaches depth z [m], and AF is the attenuation factor given by (Stenemo et al., 2007): where DT50 [d] is the half-life of the pesticide. AF expresses the fraction of the applied pesticide that will reach depth z and can take values between 0 (none of the applied pesticide will reach depth z) and 1 (all the applied pesticide will reach depth z). 355 It is well-known that organic matter concentration and microbial population density decrease with depth in soil profile, hence both the pesticide retardation and decay rates are expected to decrease with depth (Jury et al., 1987;Kookana et al., 2005). To account for this effect, the model divided the soil profile into three regions (Fig. 2): a) the A horizon (topsoil), which was assumed to be 30 cm thick and to be the most microbially active region; b) the B horizon (subsoil), which was also assumed 360 to be 30 cm thick; and c) the vadose zone, which extends from the bottom of the B horizon to the groundwater table and was assumed be microbially the least active. Furthermore, in the VDLM catchment, the presence of low-permeability soil pans was believed to be widespread due to intensive management. When such soil pan is present, it was assumed to be 10 cm thick and to be situated within the B horizon (Fig. 2). Each of the regions were characterised by uniform values of volumetric water content, soil bulk density, organic carbon content and decay rates (Appendix A). 365 where RFi and DT50i are the retardation factor and half-life in zone i, and w_i [d] is the average time it takes the infiltrating water to travel through a given zone i: 375 where w_i and di [m] are the effective volumetric water content and the thickness of horizon i, respectively.

380
The aim of the leaching risk model is to predict the pesticide load reaching the groundwater Lgw (Eq. 1). To do this, an effective attenuation factor AFgw is calculated by multiplying the AF for each of the three soil profile regions: Overall, the pesticide leaching risk assessment reflected a set of soil-specific factors (organic carbon concentration, bulk density, water content and thickness for each of the three layers), pesticide/management-specific factors (application rate, Koc and half-life), and climatic factors (rainfall and temperature from which groundwater recharge and runoff rates were estimated).

Overland flow risk assessment 390
The transport of pesticide from the site of application to surface water via overland flow is complex. Pesticides can leave a field either dissolved in runoff water or attached to eroded soil colloids. Because the amount of eroded soil lost from a field is usually small compared with the runoff volume, losses via runoff are generally considered more important than losses via erosion for most pesticides. Only for strongly adsorbing pesticides, the erosion pathway becomes the more dominant (Reichenberger et al., 2007). However, in the present study, potato cropping was the dominant land use, which has been shown 395 to be highly erosive (Vinten et al., 2014), therefore the runoff and erosion transport pathways have been lumped into one for the modelling of the overland flow losses. Other processes such as spray drift and volatilisation may also transport pesticides directly to surface water, but these were considered less important relative to the other pathways and therefore not included in the modelling.

400
The overland flow pesticide load Lof [kg ha -1 ] was calculated using Eq. 2 and assuming that the fraction of applied pesticide that will run off a given field and reach the reservoir (AFof) can be estimated as: where tro [d] is the time passing between pesticide application and the occurrence of a runoff event, and fslope, fdist and fbuffer are the slope, distance and buffer correction factors, respectively, given by: where S [degrees] is the local slope, dfr [m] is the distance from the field to the reservoir, and E [%] is a retention efficiency of 415 a field buffer strip. The above approach for modelling overland flow attenuation follows on from REXTOX (OECD, 2000).
In REXTOX, the time between pesticide application and the occurrence of a runoff event is assumed to be three days, whereas we assumed the time to depend on the month of application (see Appendix A). For simplicity and because we considered the contribution from both runoff and erosion, the available amount of pesticide available for run-off was not corrected for sorption or for plant interception. The slope correction was assumed to be a linear function up until a local slope of 10 degrees beyond 420 which no correction took place; a similar slope cut-off value is assumed in REXTOX and in e.g. Drewry et al. (2011). The buffer correction factor was informed based on a review of buffer retention efficiencies (Reichenberger et al., 2007). It should be noted that REXTOX only considers pesticide losses via runoff from fields immediately adjacent to surface waters and therefore does not include a distance correction factor. Here, we assumed that the attenuation was inversely proportional to the distance from the reservoir. The distance correction and cut-off value were based on the assumption that no significant pesticide 425 attenuation would occur within 1m of a watercourse. This assumption was informed by Reichenberger et al. (2007) who state that:' …even if pesticide-loaded runoff infiltrates into a riparian buffer strip, the groundwater table below the strip will be rather shallow (unless the stream bed is deeply cut into the floodplain), and the groundwater feeds into the nearby stream… we conclude that riparian buffer strips are most probably much less effective than edge-of-field buffer strips.'.

430
The above presented approach enabled us to a) assess relative pesticide loss risk from all fields in the study catchment, b) compare overland flow risk to groundwater leaching risk for all fields and c) evaluate optimal spatial targeting and the effect of available mitigation measures in the whole study catchment.

Model implementation and testing 435
The hybrid BBN model was constructed in GeNIe 3.0 (www.bayesfusion.com). Prior probabilities for network variables were calculated from data described in section 2.2 and Appendix A. Discrete variables were assigned a number of mutually exclusive 'states' with conditional probabilities captured in Conditional Probability Tables (CPTs). Prior probability distributions for continuous nodes were fitted to available data using the min, max, 5 th , 50 th and 95 th percentiles of the cumulative probability distribution (O'Hagan, 2012) in the SHELF package (Oakley, 2020) in the open source statistical modelling software R (The 440 R Project for Statistical Computing 4.0.1). In the SHELF package, several statistical distributions were fitted to these statistical moments and the distribution with the closest fit to the original percentiles was selected as a prior for modelling. In a few instances, where the best fitting distribution available in SHELF was not supported in GeNIe, the next best fitting available distribution was selected. For Koc and half-life, where we had limited information to inform the prior distribution, we assumed that these variables belonged to a truncated normal distribution (truncated at zero) and we considered the mean to closely 445 represent the median (i.e. the 50 th percentile). Further, we assumed that min closely approximated the 1 st percentile and the maximum the 99 th percentile and have used them as such. We rounded down/up min and max to the nearest decimal point to make them just a fraction smaller/larger than the 1 st and 95 th percentiles in the model fitting.
A discretised version of the model was then exported to R and applied at field level, using the package bnspatial (Masante, 450 2017). A discretisation method selected for each node is described in Appendix A. Discretisation was based on a mix of expert opinion (e.g., soil organic carbon and hydraulic conductivity, as well as groundwater pesticide load and the final surface water risk, which were discretised considering the likelihood of exceeding the drinking water standard concentration of 0.1 µg L -1 ), accepted values in literature (e.g., pesticide Koc and half-life), uniform cases (e.g., rainfall, temperature, and PET), and uniform interval width (e.g. depth to groundwater, distance, slope). Child nodes such as AET, infiltration rate and overland flow 455 attenuation were discretised using interpolation to ensure that conditional probabilities for the combination of parent node states (low/low, medium/medium, high/high) were meaningful. Application rate discretisation was based on equal counts but adjusted to ensure that change in application rates would result in a shift between risk classes, with the number of states maximised to allow sensitivity to change, while considering model run time. This was done empirically by progressively increasing the number of discretisation intervals and adjusting discretisation boundaries whilst checking that probabilities of 460 pesticide application rates conditional on Application Change Rate resulted in a probability density shift by one state.
Available spatial GIS layers were used as 'hard' evidence to set states for relevant nodes in the discretised model and produce spatially explicit simulations of probabilistic outcomes. Essentially, the discretised model was applied field by field using field-specific inputs from the GIS layers as hard evidence. 465 Uncertainty in the simulated outcomes in the spatial implementation of the model was evaluated by calculating the Shannon entropy index of the target nodes. The entropy H(X) for node X is defined as: ( ) = − ∑ ( ) log 2 (P(X)) (15) 470 The entropy quantifies the information content within a node and equals 0 if X is known with certainty and is maximised when is unknown (i.e. is given by a uniform distribution).
Sensitivity analysis of the discretised model was undertaken in GeNIe using the algorithm of Kjaerulff and van der Gaag (2000) 475 that calculates a complete set of derivatives of the posterior probability distributions over the target nodes over each of the numerical parameters of the Bayesian network, using the two modelled risk pathways and the combined risk as target nodes.
Euclidean distance measure, which quantifies the distance between the various conditional probability distributions over the child node, conditional on the states of the parent node, was used to calculate the strength of influence between variables (Koiter, 2006). 480 For model validation, simulated values of surface water risk (i.e., surface water pesticide concentration in µg L -1 ) in the hybrid model were compared to the limited available water quality observations for four active ingredients from month January -March (see Section 3.4). This was done by generating 10,000 random simulations for each pesticide in the hybrid BBN using forward stochastic simulation, with the random samples being drawn from the nodes' prior probability distributions. It was 485 found that 10,000 simulations were adequate for achieving stable model convergence. Note that in the hybrid model the prior probability distributions for the different nodes reflect the catchment as a whole (e.g., the prior probability distribution for the crop type node reflects the actual distribution of crop types in the catchment). This means that although spatial patterns and potential spatial correlation between inputs are not accounted for explicitly in the hybrid model (unlike for the spatial application of the model using bnspatial, where the model is applied on a field-by-field basis), the simulated surface water 490 concentrations in the hybrid model can be considered representative of average catchment conditions. Model credibility was furthermore evaluated using stakeholder feedback.

Simulated scenarios
A questionnaire was used to elicit stakeholder feedback regarding potential alterations to the management of crops/pesticides 495 and mitigation strategies from steering group members representing a grower, a regulator and a drinking water supplier in the study catchment to develop plausible pesticide mitigation scenarios. The agreed scenarios included: • Baseline risk for five active ingredients • Delayed pesticide application by one month (January to February and February to March) • Reduced application rate by 10%, 25% and 50% 500 • Additional buffering of fields to reduce overland pesticide runoff • Presence/absence of soil pan • Maximum intervention (combination of 50% reduction in pesticide application rate, management of plough pan, delayed application timing and field buffer installation) 505

Results and Discussion
A number of detailed mechanistic models have successfully simulated pesticide dynamics at plot and catchment scale (Köhne et al., 2009;Villamizar et al., 2020). However, detailed observational data required for the calibration and validation of complex models is not widely available to managers in many drinking water catchments. In this case study, in addition to sparse water quality observations, process-based modelling was hindered by the complex water transfers and limited gauging 510 at the reservoir outlet, which prevented the calculation of the hydrological balance due to lack of data on water transfers from neighbouring catchments and the productivity of the desalination plant. Difficulty with closing a water balance without considerable uncertainty is a known problem in many catchments that affects practical application of many modelling approaches (Beven et al., 2019). Therefore, to support decision making, we developed a probabilistic model of intrinsic vulnerability to pesticide pollution within a Bayesian framework to allow the assessment of intrinsic pesticide risk and inform 515 management. As pesticide risk assessment is inherently uncertain, due to many complex and poorly characterised processes, the graphical BBN model helps to improve the transparency of the risk management process (Carriger and Newman, 2012).
Furthermore, the probabilistic assessment provided by the BBN methodology is more in line with the classical definitions of risk than the more commonly used single-value risk quotients .The model represents key processes to capture combined uncertainties stemming both from observational data and limited knowledge (Sahlin et al., 2021). 520 3.1 Can we characterise the spatial and temporal variability of pesticide pollution risk from groundwater leaching and overland flow using limited observational data?
The causal structure of the hybrid BBN model designed in GeNie is shown in Figure 3. The network consists of 45 nodes and 75 arcs. The results of the spatial simulation of the groundwater leaching pesticide flux, the overland flow pesticide flux, and 525 the overall surface water risk are shown in Figures 4-6. The spatial application of the discretised model for five active ingredients has mostly shown a uniform low degree of pesticide leaching to groundwater across the 3.1 km 2 study catchment, with the exception of prosulfocarb applied at the highest application rate (Fig. 4). The largely low groundwater leaching risk is not surprising, due to the very high groundwater attenuation rates resulting from the considered pesticides being neither particularly mobile (all have relatively high Koc values) nor persistent (all have relatively short half-lives) ( Table 2). Piffady 530 et al. (2020) also found that sub-surface vulnerability was the least discriminating spatial layer, as compared to other hydrological pathways.
Conversely, the overland flow pesticide loads showed a distinct spatial variability, with most risky fields located on the steepest fields closest to surface water bodies (Fig. 5). Figure 5 also suggests that more of the risky fields are located around the west 535 stream, which agrees with the fact that the observed pesticide concentration levels in the west stream are generally greater than in the east stream (cf. Table 1). This can be explained by more fields being treated with pesticides in the western part of the catchment (i.e. more fields where the dominant land use is potato) but also by more permeable hydrogeological formations (i.e. blown sand) and soils located in the south east part of the catchment, and hence less runoff expected to be generated in this area (Fig. 1b). As the overland pesticide load is more closely related to the final surface water risk than is the groundwater 540 pesticide load (Fig. 7, Table 3), the resulting risk maps for overland load and surface water risk load look similar (Figs. 5-6).
Entropy calculation for overland flow pesticide loads (Fig. 5) and surface water risk (Fig. 6) suggests that in the baseline scenario, the risk assessment status class is more certain for the fields closest to the streams, while for the maximum intervention scenario, the uncertainty is more evenly distributed across the catchment. The assessment of groundwater leaching 545 pesticide loads is generally more certain (Fig. 4). Regardless of the absolute values, the relative difference in entropy between different management scenarios and risk status classes is informative for informing management interventions and safeguarding managers from putting too much or too little confidence in the final risk assessment (Sahlin et al., 2021).
The application of continuous and hybrid networks in environmental risk assessment is rare (Kaikkonen et al., 2020). The 550 integration of BBNs with GIS for spatial risk assessment has recently been increasing but is still limited (Carriger et al., 2021;Guo et al., 2020;Kaikkonen et al., 2020;Pagano et al., 2018). A major advantage of the BBN approach presented here over existing index methods for pesticide risk assessment is that the risk and the associated uncertainty can be determined and mapped, thereby allowing the confidence in the results to be directly assessed. The hybrid network allows more detailed characterisation of multiple processes and their uncertainty, than a typical index-based GIS method. However, the need to 555 discretise the network for spatial application currently presents a major methodological limitation, leading to loss of information, and would merit further research and development.

Which factors are most influential on intrinsic pesticide pollution risk?
Sensitivity analysis has identified the most influential parameters affecting the pesticide leaching and the overland flow pesticide loads. Figure 7 shows the result of sensitivity analysis graphically, with nodes coloured in red being more important for the calculation of the posterior probability of the pesticide risk nodes. Pesticide pollution risk was particularly sensitive to 580 crop type, time of application, overland attenuation, slope and proximity to the surface water body. Crop type directly determines the expected amount of pesticide applied in a given field, whereas the time of application affects expected rainfall and length of a dry spell following application. Alongside the soil hydraulic conductivity, the latter two influential variables determine the amount of infiltration and overland flow. It is apparent that not all variables in the model contributed strongly to the final risk assessment, hence model simplification may be possible. However, it may be advisable to test model 585 transferability to other locations first to confirm these relationships, before omitting potentially uninfluential variables. For example, depth to groundwater appears to be uninfluential in this study catchment, which may be explained either by the relatively uniform shallow depths and uncertain hydrogeological data or by most pesticide retention and degradation taking place in the A and B soil horizons, with limited influence of the vadose zone. These hypotheses could be tested in a study catchment with a better understanding of the sub-surface. Evapotranspiration calculations also appear to have limited impact 590 and could potentially be omitted from the model for simplification. While the above influential variables have been previously identified among important factors for rapid overland flow risk (Bereswill et al., 2014;Tang et al., 2012), our modelling approach allows to distinguish between generic risk factors and variables with the greatest influence on pesticide pollution risk in a local catchment setting. Coupled with the ability to identify critical source areas, and how they change at a monthly time step, our approach offers a major advancement as compared to static GIS-based risk index assessment approaches (e.g. Quaglia 595 et al., 2019). Figure 7 shows the results of the strength of influence analysis, where the thickness of the arrows represents the strength of influence between two directly connected nodes based on the Euclidean distance between the probability distributions. The top 20 most closely related variables with Euclidean distance > 0.5 are presented in Table 3. All the relationships are intuitive 600 and build confidence in reliable specification of the conditional probabilities and hence in model simulations.

What is the effectiveness of available management interventions on pesticide risk reduction?
The BBN model was applied to evaluate the effectiveness of the following mitigation measures on reducing the pesticide risks: delayed timing of pesticide application; 10%, 25% and 50% reduction in application rate; additional field buffers; and 615 presence/absence of soil pan. Figures 8-9 show the results of the simulated management scenarios on, respectively, the overland flow pesticide load of metobromuron and the groudwater leaching pesticide load of prosulfocarb (similar results for the other pesticides can be found in the Appendix B).
The time between pesticide application and the first runoff event is often considered critical for the mobilisation and loss of 620 pesticide via runoff, and hence, avoiding application in months with a higher probability of runoff events can potentially lead to a reduction in risk (Reichenberger et al., 2007). The figures suggest that delaying the application of pesticide until March results in a decrease in runoff risk for metobromuron but has limited impact on the leaching of prosulfocarb to groundwater, as groundwater risk is not related to the length of a dry spell and associated pesticide degradation following pesticide application. 625 Reduction in application rates unsurprisingly results in a reduced risk, particularly in the groundwater leaching risk (Figure 9), which is reduced to Low even after 10% reduction. This is in line with findings of Reichengerger et al. (2007) who concluded that application rate reduction was one of a few mitigation measures that could address pesticide leaching to groundwater with some confidence. However, while reduction in application rates is an easy measure to implement and may lead to cost savings, 630 it can only be implemented to a point whereby it remains effective, thus potentially limiting its acceptability as a mitigation measure to farmers (Bereswill et al., 2014). Introduction of buffers reduces the runoff risk to a similar extent as a 50% reduction of pesticide application rates for overland flow risk. Buffers, particularly edge of field buffers, have been found to be effective mitigation measures for the overland flow pathway (Bereswill et al., 2014;Reichenberger et al., 2007). Although buffer effectiveness is variable, most studies report efficiencies over 60% (Bereswill et al., 2014).In addition, permanent 635 buffers are seen as an acceptable mitigation intervention, as long as farmers can be compensated for potential losses of growing area (Bereswill et al., 2014). Managing and removing potential plough pans increases the amount of infiltration into soils, thus reducing the pesticide runoff risk to an extent that is comparable to 10% reduction in application rates (Fig. 8). By combining all available maximum interventions of 50% reduction in pesticide application rate, management of plough pan, delayed application timing and installing additional field buffers, the probability of all types of risk is notably reduced (8)(9). 640 Whilst some of these measures may be deemeed less acceptable or difficult to implement (e.g. delayed application timing (Bereswill et al., 2014)), an ambitious mitigation approach may be required to achieve real water quality improvements (Villamizar et al., 2020).

Model validation
We constructed a causal network, where the model structure was informed by expert knowledge. However, Bayesian 660 Networks can also be used as machine-learning associative tools that are suitable for deriving patterns in datasets without a specific response variable. It could be argued that pesticide risk, expressed as load or concentration of pesticide in different potential loss pathways (overland flow or groundwater leaching) is a latent variable without available observational data (Piffady et al., 2020), making it difficult to calibrate or validate a risk model. Hence, model credibility and salience (Cash et al., 2005) need to be evaluated by experts and stakeholders. Here, we implemented a simple validation approach to confirm 665 that the model predictions fall within the realms of credibility, using the limited observational data to validate the expertbased model. Figure 10 shows a comparison between the probability density distributions based on 10,000 surface water risk simulations for each active ingredient in the hybrid network and the limited observational data (in µg L -1 ) available for the months January -March between 2016 and 2019. The model typically over-estimates the simulated risk for glyphosate and pendimethalin, 670 albeit with low probability of high values. The simulated and observed distributions for prosulforcarb are comparable, whilst the model seems to under-estimate the risk from metobromuron. However, it has to be noted that the very few observations available for metobromuron (N=8) seem to be higher and less accurate than for the other pesticides. It should also be noted that the developed model was never intended to represent the complex transport and fate processes in the catchment in detail or to accurately simulate the pesticide concentration levels in the reservoir, so the comparison in Figure 10 was mainly carried 675 out as a sense check of the model predictions. Overall, model simulations appear conservative, which is helpful in terms of informing a precautionary management approach. Further model refinement could focus on constraining the upper simulation values throughout the model. A qualitative 'reasonable fit' visual inspection has been shown to be an effective means of assessing model performance using diverse incomplete data sets (Ghahramani et al., 2020). Ghahramani et al. (2020) found that the ranking of confidence in model predictions between determinands was related to data availability as much as to the 680 model itself, with pesticide simulations performing less well than those for hydrology, sediments and phosphorus.
Further validation approaches could be explored in future implementations. The spatial application in the R package bnspatial allows to simulate expected quantities, based on the median value of each discretisation interval. Hence, by multiplying the expected loads from each field with the probability of the field falling into each discretised interval, then summing the 685 resulting pesticide masses over all fields in the catchment and dividing by the reservoir volume, a concentration in the reservoir water could be estimated for each pesticide and month. This could then be compared to measured concentrations if available for further model validation. However, this deterministic calculation would be heavily reliant on the discretisation of the target node in question and, coupled with rare extreme high values generated by the stochastic model, would make the validation uncertain. Hence, it would be best applied in combination with dynamic discretisation, if available, and with further model 690 development constraining upper simulated values.

Limitations and outlook
A unique advantage of a BBN is the ability to inform probabilistic decisions on the basis of incomplete data (Panidhapu et al., 2020) and address 'what-if' counterfactual scenarios (Gibert et al., 2018;Moe et al., 2021) as well as the ability to integrate 700 data of different quality from different sources and disciplines. In machine learning, the selected mathematical approach needs to be based on the target question to be answered and be aligned with the properties of the available data (Gibert et al., 2018).
This study presents a novel approach to pesticide risk analysis that matches the question in hand with sparce data in a poorly monitored drinking water catchment.
The BBN model could easily be extended to consider other pesticides by changing the pesticide-specific properties, while 705 greater structural changes would be needed to simulate the cumulative risk from total pesticide concentrations. However, the developed BBN also has several limitations and there are various input parameters and elements that could be refined and improved. The BBN focuses only on aquatic risks resulting from intentional application of pesticide in agriculture and does not consider potential point sources of pesticide contamination (such as misuse, accidental spillages, disposal of pesticides or cleaning of application equipment). Although quite detailed in process representation that is based on established mechanistic 710 approaches, the modelling of pesticide leaching and runoff in the BBN is still simplified and could be extended to e.g., consider preferential flow pathways and/or capture the slower rate of groundwater pesticide leaching as compared to fluxes via surface runoff. Soil water balance and hydrology are critical for pesticide risk assessments, the details of which are challenging to capture with a BBN or an index-based model. Hence, future developments of the approach could include development of an improved probabilistic soil hydrology model linked to the pesticide risk model. 715 Discretisation is recognised as a major limitation of BBNs (Nojavan A. et al., 2017). Whilst here we constructed a hybrid network that has largely allowed us to avoid the loss of information associated with discretisation, this advantage was lost in the spatial application where existing mathematical and software limitations prevent direct coupling of a hybrid network with GIS. We suggest that this limitation could be an interesting and fruitful avenue for further research and methodological 720 development, e.g. by developing software applications that allow automated dynamic discretisation (Fenton and Neil, 2013) coupled with GIS. Finally, model validation would be helped by confronting with field data in a highly monitored experimental catchment, thus also allowing to evaluate the model transferability.
Notwithstanding the above limitations, this modelling approach satisfies many of the requirements of an 'ideal' model to 725 support environmental decision making, as set out by Schuwirth et al. (2019), in that it 'can be directly linked to management objectives, predicts effects of management alternatives without bias, includes adequate precision and a correct estimate of prediction uncertainty (...) and is easy to understand'. We also developed the model with the final criterion of 'easy transferability in space and time' in mind, and this could be tested in future applications.

Conclusions
In this study we present a spatial Bayesian Belief Network (BBN) that simulates inherent pesticide risk to groundwater and surface water quality, identifies critical source areas and informs field-level pesticide mitigation strategies in a small drinking water catchment with limited observational data. The BBN accounted for the spatial heterogeneity of surface water risk from pesticides, taking into account the spatial distribution of soil properties (texture, organic matter content, hydrological 735 properties), topographic connectivity (slope, distance to surface water/depth to groundwater) and agronomic practices; temporal variability of climatic and hydrological processes (temperature, rainfall, evapotranspiration, overland and subsurface flow) as well as uncertainties related to pesticide properties and the effectiveness of management interventions. The risk of pesticide loss via overland flow and leaching to groundwater were simulated for five active ingredients. Overland pesticide pollution risk from overland flow showed clear spatial variability across the study catchment, while groundwater leaching risk 740 was more uniform. The effectiveness of mitigation measures such as delayed timing of pesticide application, reduction in application rates, installation of additional field buffers; and management of soil plough pan on risk reduction were evaluated.
Combined interventions of 50% reduced pesticide application rate, management of plough pan, delayed application timing and field buffer installation greatly reduced the probability of high-risk from overland flow. The advantages of the presented BBN approach over traditional index-based methods include its ability to integrate diverse data sources (both qualitative and 745 quantitative) for a field-scale assessment of critical source areas of pesticide pollution in a data sparce catchment, with explicit representation of uncertainties. The graphical nature of the decision support tool facilitates interactive model development and evaluation with stakeholders to build model credibility; while its flexible and dynamic nature allows for performing both predictive and diagnostic reasoning based on observations, which can be linked to spatially explicit data, thus improving pesticide risk management. 750

Appendix A Model description
where d is the total soil profile depth of the topsoil and subsoil combined (assumed to be 60 cm), di and Kunsat_i are the thickness and unsaturated hydraulic conductivity of horizon i, respectively. The calculation can take presence of a soil pan into account. Soil hydrological unit and horizon specific unsaturated hydraulic conductivities are derived from the HYPRES database and depends on whether a soil pan is present (SP) or not (No SP). The hydraulic conductivity is assumed to follow a log-normal distribution with mean (µ) and standard deviation ( The time between pesticide application and the first runoff event is determined by the likely length of a dry spell. Probability of dry spell length is calculated from daily rainfall data from 2014-2019 using the method by Hills and Morgan (1981) and assuming days with less than 0.25 mm rainfall are dry. Retardation factor for the topsoil, subsoil and parent material/vadose zone. The retardation factor describes the velocity of the solute pesticide relative to the infiltrating water. Hence, a RF=1 corresponds to a solute not experiencing any retardation due to adsorption (e.g. a tracer), whereas a RF=4 means that the solute travels 4 times slower than the infiltrating water etc. For a given pesticide and soil horizon, the retardation factor is calculated as: = 1 + * * Discretisation boundaries were based on expert opinion. The attenuation factor during vertical solute transport through the topsoil (AF A). The calculation assumes 1D plug flow transport with the infiltrating water, linear retardation and first-order decay (Stenemo et al., 2007): where dA is the thickness of the topsoil, assumed to be 0.3 m. AF_A is the fraction of the applied pesticide that will reach the bottom of the topsoil and can take values between 0 (none of the applied pesticide will reach the bottom of the horizon) and 1 (all the applied pesticide will pass through the horizon). Discretisation boundaries were based on expert opinion on a logarithmic scale to reflect the skewed distribution. where dB is the thickness of the subsoil, assumed to be 0.3 m. It is assumed that the half-life during transport through the B horizon is 4 times longer than in the topsoil. Discretisation boundaries were based on expert opinion; a logarithmic scale was used to reflect the skewed distribution. The thickness of the vadose zone is given by the 'Depth to groundwater' node (Depth) minus the thickness of the topsoil and subsoil. It is assumed that the half-life during transport through the vadose zone is 1000 times longer than in the A horizon. Discretisation boundaries were based on expert opinion; a logarithmic scale was used to reflect the skewed distribution. The combined attenuation factor for the soil horizons and the vadose zone describes the fraction of the pesticide applied at the surface that will eventually reach the groundwater: GW_ = _ * _ * _ Discretisation boundaries were based on expert opinion; a logarithmic scale was used to reflect the skewed distribution. The pesticide amount leaching to groundwater is calculated from the groundwater attenuation factor (GW_AF) and the amount of pesticide applied to land (AR): ℎ = * * ℎ Leaching to groundwater is considered high if the pesticide load to groundwater exceeds 0.0001 kg ha -1 . If a load of 0.0001 kg ha -1 is mixed in the top 0.1 m of the groundwater, this will result in a concentration of 0.1 µg L -1 , which is the drinking water standard.
Medium 1.0E-5 -0.0001 High 0.0001 -7 Fraction remaining (fdecay) Very high < 0.1 Fraction of pesticide that will remain following application and decay during a dry period before the first rainfall event.
= exp (− * ln (2)  50  ) where tro [days] is the dry spell length. Discretisation was based on interpolation to ensure that conditional probabilities for combination of parent node states are meaningful. The pesticide amount reaching the reservoir with runoff is assumed to be a function of the amount of pesticide applied to land (AR), the overland attenuation factor (AFof) and the fraction of overland flow to maintain pesticide mass balance: = * *(1-fleach) The discretization is based on the same consideration as for the 'Groundwater load' node Medium 1E-5 -0.0001 High >0.0001 Surface water risk [µg L -1 ] Low <0.01 Surface water risk is the sum of the groundwater and overland flow loads. Discretisation was based on predicted likely pesticide concentration in the reservoir. This was calculated by multiplying the combined loads by the total field area in the catchment (Ac=192 ha) and dividing by the water volume in the reservoir (Vres = 938,700 m 3 ): = ( + ) * / *10 6 The risk was considered high if the resulting concentration was likely to exceed the drinking water standard for pesticides 0.1 µg L -1 .

Code and data availability
The code and data cannot not be made publicly available due to funder restrictions and privacy concerns. However, it may be available on request by contacting Director, Operations, Jersey Water, Mulcaster House, Westmount Road, St Helier, JE1 775 1DG, Jersey, UK.

Author contributions
MG led conceptualisation, funding acquisition and project administration; MT, AL, ZG and AV undertook formal data analysis; MG and MT led model development and manuscript preparation; ZG led data curation and visualisation. All authors contributed to the methodological development, manuscript review & editing. 780