Chemical exposure and risk in river basins: Bayesian integration of environmental model simulations and monitoring data

1 Introduction

1.1 Background

The Green Deal lies at the core of the European Commission's strategy to achieve a range of Sustainable Development Goals. There is growing regulatory recognition that exposure of ecosystems to mixtures (co-exposure) of chemicals may pose risks not identified through assessments of individual substances (as assessed in regulatory traditional approaches). This is intended to be addressed by the EU's Chemicals Strategy for Sustainability (https://environment.ec.europa.eu/strategy/chemicals-strategy_en) which aims to incorporate environmental risk of “combined exposure,” stemming from unintentional chemical mixtures in the freshwater ecosystems, into regulatory risk assessments (Backhaus et al., 2025).

A previous EU funded project “SOLUTIONS” (https://www.solutions-project.eu/) developed advanced and efficient tools to help identify substances and chemical mixtures that are drivers of risk. This was done through various monitoring methods and in combination of different process-based (exposure) models. Furthermore, comprehensive models and databases were integrated to enhance the understanding of chemical exposure and effects. Moreover, SOLUTIONS established a conceptual framework for the assessment, prioritization, and mitigation of pollutants and their mixtures to protect and reduce risks to ecosystems and human health. More information on the developed tools, methodology and outcomes can be found in Posthuma et al. (2019) and van Gils et al. (2020).

Within the ENCORE project (ENvironmental CO-exposure and Risk Estimation), the authors are developing a probabilistic modeling framework for the assessment of chemical risk to aquatic ecosystems at the watershed level. Chemical risk assessment can incorporate many types of information including chemical use, emission, transport, degradation, transport and fate in the environment, as well as chemical hazard (toxicity). Process-based models such as those developed in SOLUTIONs can make efficient use of such diverse information types, for simulation of chemical exposure data. However, even state-of-the-art simulation models are constrained by assumptions and quality of input data, which limits their performance under natural and variable environmental conditions. Nevertheless, environmental management of chemicals requires decisions based on the best available knowledge and assessments. Therefore, the ENCORE project is developing a framework to allow for better integration of multiple sources of information, including process-based model simulations (exposure predictions), chemical monitoring data, toxicity metrics and other types of information. This integration framework makes use of probabilistic modeling, to account for variability in the exposure information throughout the workflow. Furthermore, the probabilistic framework will provide estimates of environmental risk in terms of both probability and severity.

Environmental monitoring data (observations) can provide a snapshot of local chemical exposure to support retrospective risk assessment. However, environmental monitoring is resource-demanding and therefore limited in scope: typically covering only certain substances, specific locations, and the conditions at the time of sampling, often for a specific group of chemicals deemed important for exposure investigation (Figure 1A). Besides, the resulting chemical monitoring data can contain a considerable proportion of values below the limit of detection or quantification, which may still contribute to the total chemical risk. In contrast, prospective assessments of chemical mixtures typically rely on process-based simulation models such as environmental fate and transport models (Mentzel et al., 2022b; Oldenkamp et al., 2018; Piffady et al., 2020; Troldborg et al., 2022; Figure 1A). However, both epistemic uncertainty (stemming from incomplete or missing data) and aleatory uncertainty (reflecting natural variability) are unavoidable in exposure modeling. Moreover, exposure models might ignore legacy chemicals, which are no longer produced or used due to regulatory decisions, but which can still occur in watersheds and thereby continue to pose a risk to aquatic ecosystems (Rasmussen et al., 2015).

1.2 Probabilistic approach and Bayesian networks

Traditional procedures for environmental risk assessment of chemicals involves the calculation of scores such as a single-value risk quotient: the ratio of the predicted or measured chemical concentration to an environmental threshold value such as a hazard concentration (De Bruijn et al., 2002; Mentzel et al., 2022b; Verdonck et al., 2003). While challenges related to uncertainty are recognized in these procedures, uncertainty is usually accounted for in a pragmatic manner, for example selecting the lower 5% of a species sensitivity distribution as a toxicity metric (HC5). The HC5 is a common toxicity benchmark used in environmental risk assessment and is derived by a species sensitivity distribution that is fitted to toxicity (effect) data. It is assumed that 95% of species are protected, as the HC5 provides a conservative estimate that accounts for more variability in sensitivity across taxa than single species test estimates (De Bruijn et al., 2002). Furthermore, assessment factors are commonly applied to account for uncertainty by lowering the environmental threshold value and thereby making the risk assessment more protective. However, some of this uncertainty can be more explicitly addressed using probabilistic approaches, where risk is defined and calculated as the probability of a specified outcome (Mentzel et al., 2022a; Moe et al., 2022). Probabilistic approaches can represent all model variables, functions and parameters by probability distributions, and thereby enable the characterization and propagation of uncertainty throughout the entire model. As risk assessment becomes increasingly data-driven and AI-supported, probabilistic approaches are expected to play a growing role (Maertens et al., 2024).

Advanced statistical modeling such as hierarchical Bayesian models are becoming more common for exposure assessment in ecotoxicology [e.g., Wolf and Tollefsen (2021)]. However, monitoring data might also be used more efficiently for refinement and enhanced accuracy of predicted environmental concentrations by use of Bayesian updating of probabilities. Although Bayes' rule (Figure 1B) dates back to 1763 (Kjærulff and Madsen, 2013), Bayesian updating with new evidence is not commonly used for large-scale chemical exposure and risk assessment.

Bayesian networks (BN) are a probabilistic method increasingly used to support environmental risk assessment (Moe et al., 2021). Recently, BN modeling has also been explored for co-exposure risk characterization (e.g., Welch et al., 2024). A BN is a graphical probabilistic model that represents a set of variables and their conditional dependencies via a directed acyclic graph. The edges (links) can have a causal interpretation and a corresponding set of conditional probability distributions, which factorize the joint probability distribution over the variables of the model (Kjærulff and Madsen, 2013). In Bayesian inference, the posterior probability (prediction) is proportional to the product of the prior probability (assumptions) and the likelihood function (new evidence). In other words, the updated (posterior) BN model prediction can be considered as a compromise between the assumptions in the model and new evidence or observations (Figure 1B). A major advantage of this approach is that BNs can be constructed from tentative information sources (prior probability distributions) and updated/refined with independent data (e.g., monitoring data). Mentzel et al. (2022b) integrated predictions from a process-based pesticide exposure model into the exposure assessment and combined these with toxicity test data in the effect assessment within a BN to predict fully probabilistic pesticide risks to freshwater aquatic ecosystems. Building on this approach, Mentzel et al. (2023) as well as Martínez-Megías et al. (2023) used BNs as a metamodel to probabilistically link multiple information sources, including outputs from a case-based effect model and a process-based pesticide exposure model. Their work predicted pesticide risks to aquatic communities under varying climate and application scenarios. In another study, Oldenkamp et al. (2024) developed a BN model to represent key properties of a process-based pesticide exposure model (WISPE) in combination with a precipitation index. This BN served as an emulation model for WISPE, to be used in combination with climate projections while avoiding the computational demands of the process-based model. Collectively, these studies demonstrate the feasibility and potential of Bayesian networks for estimating pesticide risks in various freshwater ecosystems.

1.3 Objectives

Within the ENCORE project we are developing a probabilistic modeling framework which comprises two main components: a process-based exposure model for simulation environmental concentrations, and a Bayesian network (BN) for predicting exposure and risk, both for individual substances and mixtures. The goal of this probabilistic framework is to identify chemicals with high contribution to mixture risk, by accounting for spatial and temporal exposure patterns in watersheds across large regions of Europe. The exposure model (referred to as “the ENCORE exposure model” in the following) builds upon a process-based chemical exposure model originally used in the SOLUTIONS project (Lindim et al., 2016; van Gils et al., 2020). An important purpose of BN developed in the ENCORE project is to integrate simulations from the process-based model chain (representing prior probabilities) with new monitoring data (representing new evidence; Figure 1). The aim is to integrate and propagate both the key signals and uncertainties from the process-based model as well as from other information sources (Table 1), throughout the BN. Compared to the more complex process-based model (Section 2.2), the BN will have lower resolution and precision, but will be more streamlined and can be run more efficiently. More specifically, the study reported here explores the use of BN methodology to update prior probabilities with new evidence from chemical monitoring data. This exercise aims to improve the precision and accuracy of both model prediction and parameters. Efficient use of geographical information systems (GIS) tools in combination with probabilistic modeling will ensure the different data sources are effectively integrated at appropriate spatial scales. This applies to the exposure part of the risk assessment, while for the toxic effect part we currently use single values.

Table 1. Overview of main data sources used in this study.

“Exposure model” refers to the process-based simulation covering chemical use, emission, transport, fate and concentration in subcatchments (Section 2.2). WISE, Water Information System for Europe; WFD, Water Framework Directive; BN, Bayesian network; RBD, river basin district; RQ, risk quotient (Section 3.2).

The objective of this paper is to describe the development, exploration, and potential use of the BN model, covering all steps of the workflow from data compilation and processing to probabilistic risk visualization and interpretation. We have developed the workflow so far based on a pilot study in Belgium and selected a subset of 15 pesticides, with the longer-term vision of upscaling the approach to larger regions of Europe.

2 Data

2.1 Pilot study: pesticides in Belgium

The pilot study area for the BN model development is Belgium (Figure 2). The process-based exposure model applied a spatial structure with the country divided into 144 subcatchments. In the BN, the spatial extent of the pilot study was constrained to five river basin districts (RBDs), based on considerations of their data availability, with the following international codes: BEESCAUT_RW, BEMAAS_VL, BEMEUSE_RW, BERHIN_RW, BESCHELDE_VL. In terms of administration and water management, these RBDs represent the river catchments Escaut/Schelde (Scheldt), Meuse/Maas and Rhin (Rhine), in the Belgian regions of Wallonie (RW) and Flanders (VL), respectively. Since delineation of the subcatchments is not fully aligned with delineation of the RBDs, the analysis was made for subcatchments that fall within respective RBD with more than 50% of their area; if a subcatchment's area is shared between more RBDs, the whole subcatchment is assigned to the RBD that takes the largest share of respective subcatchment. Of 144 subcatchments analyzed in total, 30 have a significant share of area (more than 10%) outside the assigned RBD.

For the pilot study, 15 pesticidal substances were selected with relatively high application rates in Belgium (Table 2). Pesticides were modeled using an area-specific approach focused on agricultural land. The SOLUTIONS agricultural land use dataset was used to include spatially explicit land use (e.g., crop type coverage and seasonal use) combined with area-specific application rate maps developed by Porta et al. (2025). In that study, pesticide application rates were estimated at a spatial resolution of 250 m across the EU. The resulting application rates were subsequently linked to the subcatchments used in the spatial schematization of the ENCORE exposure model.

Table 2. Selected chemical substances with HC5 values (hazardous concentrations to 5% of species) and assigned substance groups (Pistocchi et al., 2019).

Several criteria were considered for the selection of the 15 active pesticide substances used in the pilot study. First, we relied on the harmonized, EU-wide application rate maps provided by Porta et al. (2025), which compile national sales data and crop-specific use information. Substances were therefore primarily selected from those for which robust, spatially explicit data were available for the region of interest (Belgium). Second, monitoring data availability played a key role in the selection. Only substances for which consistent monitoring data were available in the WISE-6 Waterbase database (EEA, 2024) for Belgium were retained. As a reality check, concentration predictions from the ENCORE exposure model would need to be compared to field-monitoring data for the substances of interest.

The pesticides were assigned to three different substance groups, according to their chemical nature and functional mode of action (based on Pistocchi et al., 2019): fungicides (including bactericides), herbicides (including haulm destructors and moss killers), and insecticides (including acaricides). This grouping reflects differences in application practices and intended targets, which act on specific components and developmental stages within the ecosystem.

2.2 Simulated data—ENCORE exposure model

Simulated chemical concentrations for the pilot study were obtained from the ENCORE exposure model, which is a process-based stochastic simulation model covering chemical use, emission, transport, fate and concentration in subcatchments. The ENCORE exposure model can produce exposure profiles in surface waters across Europe, for individual substances representing three main chemical categories: pesticides, pharmaceuticals, REACH-registered industrial chemicals, including both legacy pollutants and currently used chemicals. Overall, the currently available model simulations cover daily concentrations of various substances in 22 728 subcatchments across 42 countries for the period 2004–2013, each with 100 stochastic runs (Monte Carlo simulations). For this pilot study, we selected only a subset of the simulations representing the 15 substances (all pesticides; Table 2) in the 144 subcatchments of one country (Belgium), for 1 year (2012), and a single model run (deterministic).

The ENCORE exposure model incorporates several methodological refinements compared to the original model from SOLUTIONS (van Gils et al., 2020). The Monte Carlo simulation framework was developed to represent uncertainty in key input parameters, such as emission factors, and to quantify how this uncertainty propagates through the model to influence predicted concentrations and eventually associated risk estimates. The model was calibrated using European hydrological datasets to ensure realistic flow dynamics and spatially consistent transport patterns.

The chemical concentrations in surface waters were simulated through two main computational steps for each individual substance. In the first step, chemical emissions were simulated by assigning substance-specific source categories and emission pathways to each sub-catchment. Distinct emission mechanisms were defined for the different chemical categories (pesticides, pharmaceuticals and REACH substances). For pesticides, agricultural application was identified as the dominant source: products were assumed applied either to soil near plant roots, resulting in emissions via runoff and leaching, or to foliage, where a fraction is taken up by plants and the remainder is transferred to the environment. These emission pathways simulated the transfer of chemical mass across environmental compartments and calculated the amount of each substance entering surface waters within each subcatchment. In the second step, the water quality module computed chemical concentrations in the aquatic environment by explicitly accounting for hydrological connectivity between subcatchments. In-stream processes such as dilution, degradation, and phase partitioning were also represented. Together, these processes determined the spatial distribution and magnitude of chemical concentrations across the sub-catchment.

The simulated exposure concentrations obtained from the ENCORE exposure model outputs at the sub-catchment level were subsequently processed (Section 3.1) for use as prior probabilities in the BN model.

2.3 Monitoring data—WISE-6 data

This project aims to make use of monitoring data on chemical substances to leverage the exposure model simulations. As a starting point, we have prioritized chemical monitoring data from the European Environment Agency's Water Information System for Europe (Waterbase—Water quality; EEA, 2024), which contains observations of hazardous substances from more than 50 million samples across Europe. This data source will be referred to as WISE-6, which is the label of the underlying EEA data flow.¹

The WISE-6 dataset was downloaded on 29.08.2025 as a set of two SQLite databases: Waterbase_v2024_1_WISE6_DisaggregatedData.sqlite² (containing chemical concentrations for samples by monitoring site and date) and Waterbase_v2024_1_WISE6.sqlite³ (containing additional spatial information such as coordinates and water body information). The WISE-6 monitoring data has the same spatial data structure as the WISE-WFD data, i.e., data from the River Basin Management Plans 2022 reporting according to the Water Framework Directive.⁴

The initial processing of the WISE-6 data (Figure 3 step II), such as filtering by country code (BE), monitoring years (2004–2013) and determinands (Table 2), was performed with the tool SQLiteStudio (https://sqlitestudio.pl/). The purpose was to reduce the data set to a size that could be further handled by standard data processing software such MS Access and R. The substances were identified by CAS code, and all substances from the fate simulations except one (Dimethenamid(–p)) were found in the WISE-6 code list.⁵ The resulting dataset for this pilot study contained 7,669 records representing in total 7 RBDs, 84 river water bodies, 89 monitoring sites and 834 sampling dates. After filtering the selected 5 RBDs, the dataset contained 7,512 records representing in total 82 river water bodies, 87 monitoring sites and 834 sampling dates.

The spatial objects of exposure model (subcatchments, Section 2.2) were mapped to the spatial objects of the monitoring data (RBDs) by GIS. If the spatial extent of a subcatchment overlapped with more than one RBD, the subcatchment was assigned to the RBD that contained the highest proportion of the subcatchment's area.

3 Methods

3.1 Workflow and data processing

The workflow for this study's probabilistic framework to assess chemical risk within a specific geographic region, involves several key steps and integrates multiple data sources. A conceptual (qualitative) model was developed, incorporating the necessary properties relevant to the region and chemicals of concern (Figure 3).

The modeling workflow consists of several key steps (Figures 3A–E). Step A begins with the collection of the main data sources (Figure 3A), which include the following (see Table 1): simulated data from a processed-based model (ENCORE exposure model)—predicted exposure concentrations that; Environmental Threshold Value (ETV) i.e., toxic effect data, in this case HC5 values (see Table 2); information on chemical groups (see Table 2), and relevant spatial information (see Figure 4). In step B (Figure 3B), the ENCORE exposure model output is processed together with spatial and toxicity data, and their empirical relationships are quantified. For the BN development (Figure 3C), these relationships are used for the discretization of variables (grouping of continuous values of the nodes into intervals), derivation of prior (initial) probabilities and parameterization of the BN with conditional probability tables (CPT) that define the relationships between nodes (see Figure 5). The completed BN model (Figure 3D) can be run based on predefined scenarios, and the outputs are visualized and summarized to support interpretation. Following this process, specific components of the BN can be identified for further refinement based on evidence. This may involve the integration of additional data (observations), in this case monitoring data (Waterbase Water quality) to enhance the model's accuracy (Figure 3E). Bayesian updating techniques are applied to improve the resolution of the model predictions (see Figure 1B; further description of Bayesian updating can be found in Section 3.3).

The spatial resolution of the BN model developed for this pilot study (RBD) was based on spatial structure of the WISE-WFD data model (Table 1). In comparison, the simulated data from the ENCORE exposure model had a higher spatial resolution—subcatchments (Figure 4-V1), derived from the E-Hype model from SMHI (Swedish Meteorological and Hydrological Institute). The reason for selecting the WFD spatial structure for this pilot study was two-fold: to simplify the workflow, considering the integration with monitoring data (Figure 3), and to enhance the relevance of the BN model predictions to water management.

The simulated exposure data needs to be aggregated along with both the spatial and temporal dimensions to allow integration of these data with other information sources with different resolutions (Figure 1). The optimal level of aggregation should preserve the main spatio-temporal patterns of chemical co-exposure and risk, while still allowing for efficient construction and running of the BN model (described in Section 3.2). The temporal and spatial dimensions considered in this study are specified in Figure 4, with alternative versions of aggregated datasets labeled V1–V8. The temporal dimension spans from days (V1—based on simulation outputs), to years (V4), while the spatial dimension is defined through subcatchments (V1—based on simulation outputs), RBDs (based on WFD; V5), country (V6) and continent (V7). These dimensions can be combined in various ways. For the pilot study of Belgium, we used the spatio-temporal units of 5 RBDs and 12 months for aggregation (V8). This information is used during the data processing approach throughout the workflow (Figure 3): within each spatio-temporal unit, the simulated concentration values were aggregated by probability distribution. For the monitoring data, the further data processing is described in Section 3.3.

3.2 BN model construction

The BN model is designed to calculate environmental risk of an individual substance by combining chemical emission and exposure data with an environmental threshold value (Figure 5A). While the overall structure is identical for each of the 15 substances (Table 2), some details are substance-specific: the discretization of continuous variables [nodes related to emission, concentration and Risk Quotient (RQ); Supplementary material S1; Table 3] and the statistical expressions for the node for “Emission to subcatchments” (Supplementary material S2). A full model documentation is provided for one example substance, diuron (Supplementary material S3). In the following, the description of the refers to the substance-generic structure (alternatively to the example of parameterization for diuron). The current version of the BN has 17 nodes representing six main types (see Figure 5B): categorical, numeric, intervals (all discrete); function (single value), continuous and Boolean.

Table 3. Table 3 Overview of nodes in the Bayesian network (Figure 5).

A central component of this study is the node “Conc.” (chemical concentration in μg/L, 12 intervals). The parameterization (definition of probability values) of this node is currently based on the spatio-temporal patterns of the simulated concentrations. The node “Conc.” is linked to the spatial node “River Basin District” (5 states) and the temporal node “Month” (12 states), which are combined into the node “Spatial & temporal settings” (60 states), via its parent node “Emission to subcatchment”. This parent node currently serves as a placeholder for further model refinement, and to enable exploration of emission scenarios. For simplicity, the “Emission” node is currently set identical to the “Conc.” node. The emission level node is defined with five states, ranging from 50% to 150% of a baseline emission level. This node can reflect potential variations in pesticide application or release, reflecting for example policy targets or worst-case scenarios.

For the calculation of the environmental risk, the BN combines the probability distribution of “Conc.” with the substance-specific “HC5” value, which represents the Environmental Threshold Value (ETV) for the given substance. The resulting RQ (15 intervals) is calculated as the concentration divided by the HC5. To simplify its interpretation, the RQ is also summarized into five levels with fixed boundaries, providing an overview to enable the comparison of RQ across substances. Finally, the Boolean node “RQ exceeds threshold?” (2 states) serves as a simplification of the risk characterization into a binary true/false statement. This node summarizes the probability of RQ intervals above the given “Threshold RQ” (6 states). By default, the value of “Threshold RQ” is set to 1, in line with conventional practice in risk assessment (demonstrated in Mentzel et al., 2022b). However, alternative levels of the “Threshold RQ” can be explored, for example, to compare risk of pesticides in months with high or low application.

The five continuous nodes labeled “MEC log” (Measured Environmental Concentrations), represent the monitoring data, and play a role in the Bayesian updating (Section 3.3). The software used in this study (HUGIN Expert) allows continuous nodes to be quantified by a normal (Gaussian) distribution with parameters μ (mean) and s2 (variance). However, environmental concentrations can be better characterized by a lognormal distribution. Therefore, the intermediate node “Conc. log,” which is an ln-transformed version of “Conc.,” serves as a parent node for the continuous nodes.

BN model predictions can be sensitive to the discretization of key variables, which should therefore be performed carefully. When considering the environmental risk of individual substances, the peak concentrations are of most interest. However, the lower concentrations ranges should also be taken into account in the discretization, since the substances with low concentrations could still contribute to mixture risk, if the BNs for individual substances (Figure 5) are combined for calculating mixture risk. Hence, an appropriate discretization is needed for both low and high concentrations. Therefore, for each substance, the discretization for all concentration-related nodes as well as for the RQ node were generated by an algorithm accounting for the substance-specific concentration range.

The central node “Conc.” was discretized into 12 intervals. This selected number of intervals balances the resolution of the variable with the need for information: a higher number of intervals would give more precise concentrations estimates, but would also require more information for parameterization (the term “parameterization” here refers to the quantification of a node's CPT, which corresponds to parameters of a process-based model). The discretization of “Conc.” is based on the range (Min and Max) of simulated concentrations, as follows. State 1 is set to the interval [0 – 0] μg/L. For states 3–11, to cover all simulated values with a margin, the range for discretization is set from 0.9 × Min to 1.1 × Max. This range is divided into 9 intervals that are equidistant in log (ln) scale. Finally, state 2 is set to the remaining interval [0 – 0.9 × Min], while state 12 gets the interval [1.1 × Max – Inf].

The intervals of the node “Conc. log” simply have the ln-transformed boundaries of the node “Conc.” Likewise, the intervals of “RQ” simply have the boundaries of “Conc.” divided by the HC5 value.

3.3 Bayesian updating

One purpose of the BN model presented was to enable Bayesian updating techniques for synthesizing environmental model simulations with new monitoring data (Figure 1A). In brief, the method makes use of Bayes' theorem to calculate a probability of a hypothesis, given the current evidence, and update it as more information becomes available (Figure 1B). This way, Bayesian updating can make use of new data in an efficient manner: instead of using the new data as model input and re-running the model, the new data can be used as evidence in combination with the current (prior) probability distributions (Figure 5B) to produce an updated (posterior) probability distribution. In comparison, for a process-based model, data-assimilation methods might be utilized to account for new measurements and improve the model performance, as an alternative to re-running the model.

When applied to a BN model, Bayesian updating (Figure 1B) can represent two different processes. In this project, we are exploring the use of these two types of Bayesian updating at different levels of spatial aggregation, for different purposes. (1) Inference: use of evidence for calculating updated model predictions, i.e., the probability distribution of a variable (Figure 5B). This simpler type of updating is based on probability calculations and will not affect the BN model structure or future model predictions. (2) Learning: use of evidence, usually multiple cases, for updating model parameters, i.e., the CPT of a node. This more advanced type of updating is based on machine-learning algorithms, in this case the expectation-maximization algorithm (Kjærulff and Madsen, 2013; Lauritzen, 1995). This type of updating will permanently change the model parameterization, and thereby affect any subsequent model predictions. Such machine learning procedures can also encompass learning of model structure (e.g., presence and direction of links), but this is beyond the scope of the current pilot study.

In this study we used the pesticide diuron as an example for both types of updating, based on the data availability and range of values (see Supplementary material S3). The initial filtering according to the spatial and temporal extent of the pilot study (described in Section 2.3) resulted in 2,167 observed values of diuron. This selection included 841 records flagged as below the limit of quantification (LOQ), for which the respective LOQ was used as the value (ranging from 0.0018 to 0.03 μg/L). The resulting dataset was further processed into two separate data files for the two types of Bayesian updating.

For the first purpose—inference (refined predictions at the water body level), only data from 2012 were selected, to match the temporal representativity of the BN model. The resulting 837 observations of diuron were organized into a table with 106 records representing different cases for inference (unique combinations of year, month, and water body), with five columns representing different 7-day periods within each month. These five columns correspond to the five nodes “MEC log 01”—“MEC log 05” (Figure 5A). However, each resulting case had only one observation, meaning that the 2012 monitoring scheme for diuron contained no more than one sample per month and water body.

For the second purpose—learning, all monitoring data from the full simulation period 2004–2013 and from the whole country (BE) were included to update the CPT of the BN using a set of cases. To maximize the number of cases, the cases were defined by individual monitoring sites rather than water bodies. The resulting number of diuron observations (still 2167) were similarly organized into a table with 2,084 records representing different cases (unique combinations of year, month, and monitoring site), with five columns representing temporal spread. For this subset, some of the cases had two observations per month and site, representing multiple sampling events in different weeks within a month.

4 Results (preliminary)

4.1 Predicted environmental exposure

Temporal patterns in the predicted exposure are evident in the aggregated ENCORE exposure model output. Figure 6 illustrates the temporal patterns for three selected pesticides in the RBD BERHIN_RW as an example. Over the course of the year, lower predicted concentrations are consistently observed for all selected chemicals (thiophanate-methyl, dicamba, and dimethoate). This temporal pattern is present across all RBDs, although the magnitude of variation differs among them. The largest differences between maximum and minimum exposure occur in BERHIN_RW (Figure 6), while the smallest differences are observed in BEESCAUT_RW and BEESCHELDE_VL (Supplementary Figure S1 in Supplementary material S4).

4.2 Predicted environmental risk

The plots in Figure 7 display the temporal patterns in risk predicted by the BN for each pesticide in the RBD BERHIN_RW (and likewise for the remaining RBDs in Supplementary Figure S2 in Supplementary material S4). These figure indicate that the patterns of temporal and spatial variability of the BN output reflect those in the ENCORE exposure model simulations (Figure 6). Among the different chemical groups, the primary risk drivers are chlorothalonil for fungicides, pendimethalin for herbicides, and chlorpyrifos for insecticides. In the case of chemical mixtures, chlorpyrifos constitutes the dominant risk driver, as it exhibits the highest individual risk contribution. The relative importance of these risk drivers varies both seasonally and across RBDs: peak risks are observed in February in BEESCAUT_RW, in July in BEMAAS_VL, in April in both BEMEUSE_RW and BEESCHELDE_VL (Supplementary Figure S2 in Supplementary material S4), and in August in BERHIN_RW (Figure 7). Overall, the results suggest that each season comprises periods of both elevated and reduced risk, with winter generally associated with slightly lower risk levels.

4.3 Updated predictions and parameters

The two types of Bayesian updating described in Section 3.3 were explored for different levels of spatial aggregation.

Inference: updating of BN predictions at low spatial resolution (RBD level) with observations at higher spatial resolution (water body level), for improved accuracy of exposure characterization. This type of updating is applicable for spatial objects (e.g., water bodies) where monitoring data are available.

An example of combining data sources with Bayesian inference is shown for exposure concentrations of the substance diuron (Figure 8). The selected case is the water body “ZEESCHELDE I” (waterBodyIdentifier “BEVL11_40”), located in the RBD “BESCHELDE_VL,” for July 2012. The BN (which represents year 2012) is instantiated accordingly for Month = 7 and RBD = BESCHELDE_VL (corresponding to Scenario 55). The updating was performed with log-transformed observed values [log(0.065) = −2.74] and variables (node “Conc. log”), for better approximation of the assumed lognormal distribution for concentrations. The prior distributions values are shown in Figure 8A for the nodes “Conc. log” (expected value μ = −26.89, variance s² = 2,933.57) and “MEC log 01” (μ = −7.58, s2 = 81.62). The five MEC nodes had identical prior distributions, given by their CPT Table 4. After the evidence for node “MEC log 01” was entered (Figure 8B), the posterior distribution of “Conc. log” was both more precise (variance s2 = 190.2 being reduced) and more accurate (μ = −4.17 being closer to the observed value). After the updating, the modal state of “Conc. log” shifted from interval #2 (−230 to −12.31, with 19.45% probability) to interval #8 (−2.954 to −1.082, with 67.4% probability), which also contains the observed value (−2.74).

Table 4. Conditional probability table (CPT) for the continuous nodes “MEC log 01” to “MEC log 05.”

Following the updating, the expected value of new observations (“MEC log 02”) also shifted closer toward the evidence and with lower uncertainty (μ = −3.57, s2 = 10.03; Figure 8C). If more observations with similar values are included as evidence, the precision of the predicted concentration is likely to improve further.

Learning: combination of CPTs representing low spatial resolution (RBD level) with observations at even higher spatial resolution (monitoring site level), for refined model parameterization. This type of updating can also affect spatial objects (e.g., water bodies) without monitoring data and potentially give more precise and/or accurate model predictions also for these objects.

An example of combining data sources by BN parameter learning is shown in Figure 9, still with diuron concentrations as the example. Learning was applied to the node “Emission to Subcatchments” (since the child nodes “Conc.” and “Conc. log” currently have only deterministic CPTs). The parameter learning exercise clearly demonstrates that the updated CPT (Figure 9B) has higher modal concentration interval than the original CPT (Figure 9A). For example, for parent state no. 1, the modal concentration interval after updating has increased from [0–4.5E-6] μg/L to [8.03E-3–0.05] μg/L. The change in variance is less conspicuous in Figure 9, but in general the variance of each column has decreased after the updating. For example, for state no. 1, the calculated standard deviation of the conditional distribution has decreased from 13.6 to 8.2 μg/L. This means that a BN with the updated CPT in this case will provide more precise predictions (with lower variance) of chemical exposure.

In this example, for simplicity, all records of diuron reported to be below the limit of quantification were represented by their respective LOQ value. This procedure is likely to overestimate the chemical concentrations, compared to other conventional procedures such as using LOQ/2, or more advanced procedures such as modeling the distribution of values below the LOQ. While a more precise estimation of values below LOQ is important in the context of mixture risk assessment, it is beyond the scope for this example, where the purpose was to demonstrate the two types of Bayesian updating.

5 Discussion

5.1 Evaluation of the Bayesian network approach for synthesis of environmental information

The use of BN methodology for assessment and management of contaminants in aquatic ecosystems have increased over the last decades (Moe et al., 2021). Many examples focus on nutrient pollution (phosphorus and nitrogen) from fertilizers in agricultural catchments (e.g., Barton et al., 2008; Negri et al., 2024; Radomyski and Ashauer, 2022; Wei et al., 2024). More recently, BN approaches have been applied for environmental assessment of other chemical stressors such as pesticides (e.g., Mentzel et al., 2023; Piffady et al., 2020; Troldborg et al., 2022) or pharmaceuticals (Welch et al., 2024). General considerations on the use of BN methodology for assessment of environmental pressures in river basins are discussed by e.g., Glendell et al. (2025). Landis (2021) reviewed several examples of BN models for environmental risk assessments with the Relative Risk Model (BN-RRM). In all of these examples, the interaction with stakeholders and decision makers improved the ability to gather datasets, establish endpoints, and produce results to inform the decision-making process. Moreover, the graphical nature of the BN models and the mapping of relative risk in the study areas facilitated the communication between risk assessors and stakeholders (Landis, 2021).

The pilot study reported here provides an example of how BN modeling can be used to synthesize environmental information from different sources and scientific domains. This pilot study has focused on the integration of three main types of information: toxicity metrics derived from experimental data, chemical exposure data derived from process-based model simulations, and environmental monitoring data from published European datasets. The BN construction enabled a representation of the main spatial and temporal patterns of the simulated chemical exposure profiles by statistical distributions (priors). Furthermore, the Bayesian updating methodology allowed for efficient inclusion of site-specific monitoring data (evidence), resulting in refined exposure profiles (posteriors) for the example substance diuron.

This pilot study demonstrates several strengths of using BN methodology for integrating large data sources for environmental risk assessment.

• (1) Representation of environmental states (e.g., chemical concentrations) by probability distributions (Figure 5B), which represents both the expected value and the uncertainty.
• (2) Representation of processes and causal links (e.g., from chemical use via emissions to concentrations) by conditional probability distributions. The distributions incorporated in the pilot study are so far based only on the simulated exposure data (Figure 6) across the spatial and temporal units (RBDs and months, Figure 4), with emission-related nodes set temporarily as placeholders without influencing the model. This part of the BN can easily be further developed to represent both emission processes and other relevant factors, by incorporating additional variables from the ENCORE exposure model.
• (3) The calculation of risk in terms of both probability and outcome (here: threshold exceedance), which is in line with definitions of environmental risk (Moe et al., 2022), although risk characterizations in practice are often reduced to single-value metrics.
• (4) The possibility for updating and refining predicted probability distributions with new evidence (Bayesian inference, Figure 8). This means that additional data sources can be efficiently made use of toward the end of the workflow, without the need for re-running or re-calibrating large models.
• (5) The possibility for updating and refining conditional probability distributions with new evidence (parameter learning, Figure 9). While this procedure is less straight-forward and more data-demanding, it can give wider opportunities for more precise assessments, e.g., for water bodies without relevant monitoring data.

The three first points above could in principle be covered by probabilistic output from the ENCORE exposure model, although it would be more time-consuming than with the BN, which can serve as an emulation model. The ENCORE exposure model is developed as a stochastic model, meaning that multiple Monte Carlo (MC) simulations with identical settings will generate different outcomes (different single-value concentrations), which can be summarized by a probability distribution. In contrast, the BN presented here can be defined as a deterministic model, i.e., multiple model runs with the identical setting will generate identical outcomes (identical probability distributions).

The deterministic property of the BN is enabled by the discretization of all variables except the final child nodes (Figure 7A). This property also allows for other functionalities such as diagnostic inference or backward calculation (e.g., Moe et al., 2020). For example, setting the condition that RQ should not exceed 1, for a given substance and spatio-temporal settings, the BN can provide an updated probability distribution on chemical use and emission. Another feature facilitated by this deterministic property is sensitivity analysis by value-of-information analysis (Mäntyniemi et al., 2009). This feature can prove useful in further BN development based on more information from the ENCORE exposure model, for example, for identifying the variables with the most influence on the exposure of different substances.

The discretization of variables in BNs also implies constraints, such as loss of precision for example in concentration variables (Figure 5B). For this pilot, concentrations were represented with lognormal distributions, which is a common assumption for environmental chemistry (e.g., Mentzel et al., 2022a), fitted to the simulated concentration for each substance, RBD and month (see Figure 9A). However, given the temporal variation in pesticide use, more specific expressions such as zero-inflated and bimodal distributions may be needed for better characterization of exposure profiles.

The spatial scale used in this study (RBD) was sufficient for displaying regional patterns in pesticide risk (Figure 7) but may be too coarse for more specific risk assessments. Further BN model development can therefore explore the use of higher spatial resolution, for example RBD subunits. Moreover, the assignment of the original spatial units (here: subcatchments) to the aggregated spatial units (here: RBD) could be refined: instead of assigning the subcathcment to the dominating RBD, the proportion of spatial overlap with each RBD could be taken into account in the aggregation. Nevertheless, information on chemical use and emission is often available only at lower resolution (e.g., country or RBD), which will constrain the resolution of the prior probability distributions.

The spatial extent of this pilot study (Figure 2) comprises only a small proportion of Europe, which is ultimately the spatial scope of the ENCORE project. Nevertheless, the BN structure developed for the pilot study (Figure 5) incorporates the multiple spatial and temporal scales (Figure 4), which will also be relevant for BN models scaled up to larger spatial regions. The selected substances (Table 2) represent only one of the three chemical categories to be covered by the project (pesticides, pharmaceuticals and REACH substances). However, pesticide exposure can be expected to display the highest spatio-temporal variation (Figure 6), since the chemical use patterns follow regional variation in growing season for various crops. Moreover, the selected pesticides represent substance groups with different MoA (Figure 7), which will determine their role in the forthcoming mixture risk calculations. Therefore, we conclude that the approach developed for integrating information in this pilot study (Figure 1) is a good starting point for extended BN development to represent larger spatial, temporal and chemical domains.

The BN presented here (Figure 5) represents a single substance, which can be efficiently combined with other substance-specific BNs as sub-models (classes) in a multi-level (object-oriented) BN. The nodes with gray outlines indicate the relevant nodes for connections to other sub-models at substance group or mixture level, for example “Spatial & temporal settings” (input) and “RQ” (output). Such a modular construction will enable an efficient and flexible construction of a multi-substance BN for mixture risk calculation. The node “Spatial & temporal settings” currently has 60 states representing all 5 x 12 combinations of RBDs and months. For spatial expansion to include more countries, additional spatial units can simply be added to the node “River Basin District.” However, to more efficiently represent spatio-temporal exposure patterns across larger regions of Europe, the statistical expressions (Figure 9A) should be refined for example by a mixed effect analysis.

The Bayesian updating of conditional probabilities demonstrated in this study (Figure 9B) makes use of artificial intelligence (AI) through a well-established machine-learning technique (expectation maximization; Kjærulff and Madsen, 2013; Lauritzen, 1995). For expansion of the BN structure to include more information from the process-based model (e.g., environmental variables), further use of this technique can be considered also learning of model structures (e.g., causal directions between nodes). Additional use of more novel AI methodologies can also be explored, as recommended for ecotoxicology (Maertens et al., 2024), for example for data mining where data gaps are identified. The monitoring data source used in this pilot study (EEA, 2024) contained observations for only four of the 15 substances. For an expanded study, additional sources for large-scale exposure data can for example include the curated database NORMAN EMPODAT (https://www.norman-network.com/nds/empodat/), containing emerging substances that are known to be present in the environment but not yet included in routine monitoring programmes. Moreover. independent regional or more local data sources can be useful for validation of BN predictions for new case studies, for example the Danube River Basin Database (https://wq-db.icpdr.org/) for case studies within this river basin.

5.2 Conclusions and outlook

This pilot study has demonstrated the development of a probabilistic model for integrating different data sources—model simulations, toxicity metrics and monitoring data—for large-scale environmental assessments. The development and use of a BN for this purpose serves as proof of concept, and a starting point for expansion of the approach to larger spatial regions. Moreover, the single-substance BNs presented here can serve as building blocks for a multi-level BN for the assessment of combined risk within and across substance groups.

Further development of the BN presented here should focus on (1) incorporating use and emission scenarios, (2) refining statistical distributions in exposure and (3) validation of predicted concentrations by independent data sources, potentially supported by further use of AI methodology, (4) incorporate probability distribution for effect assessment (see Mentzel et al., 2022b, as an example for fully probabilistic risk assessment). Further model development and interpretation will also be guided by feedback from the stakeholders involved in this project.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffwsc.2026.1789703/full#supplementary-material

Pilot 2 substance discretization. Supplementary material at the source

Pilot 2 substance priors. Supplementary material at the source

Network documentation substance diuron. Supplementary material at the source

Temporal and spatial variation for all RBDs supplementing Figures 6 and 7. Supplementary material at the source

Footnotes

Edited by: Daniel M. Perkins, Brunel University London, United Kingdom

Reviewed by: Asude Hanedar, Namik Kemal University, Türkiye

Janne Juntunen, Finnish Environment Institute, Finland

Maria Kalyva, Vitenskapskomiteen for mat og miljø, VKM, Norway

Footnotes

https://cdr.eionet.europa.eu/help/WISE_SoE/wise6

https://sdi.eea.europa.eu/datashare/s/3JiTia3qePyGxyA/download?path=%2F&files=Waterbase_v2024_1_WISE6_DisaggregatedData.sqlite

https://sdi.eea.europa.eu/datashare/s/3JiTia3qePyGxyA/download?path=%2F&files=Waterbase_v2024_1_WISE6.sqlite

https://cdr.eionet.europa.eu/help/WFD/WFD_780_2022

https://cdr.eionet.europa.eu/help/WISE_SoE/wise6/WISE6_ObservedProperty_QC_reference.xlsx