Using weighted expert judgement and nonlinear data analysis to improve Bayesian belief network models for riverine ecosystem services

(capturing con ﬁ denceamongex-pert groups), legislation and published literature. The BBN was calibrated with analyses of national monitoring data (including non-linear relationships and ecologically meaningful breakpoints) and expert judgement. We used a novel expected index of desirability to quantify the model outputs. Secondly, we applied the BBN to three case study catchments in Ireland to demonstrate the implications of changesin stressorlevelsfor ecosystem services in different settings. Four out of the seven signi ﬁ cant relationships in data analyses were non-linear, highlighting that non-linearity is common in ecosystems

Rivers are a key part of the hydrological cycle and a vital conduit of water resources, but are under increasing threat from anthropogenic pressures. Linking pressures with ecosystem services is challenging because the processes interconnecting the physico-chemical, biological and socio-economic elements are usually captured using heterogenous methods. Our objectives were, firstly, to advance an existing proof-of-principle Bayesian belief network (BBN) model for integration of ecosystem services considerations into river management. We causally linked catchment stressors with ecosystem services using weighted evidence from an expert workshop (capturing confidence among expert groups), legislation and published literature. The BBN was calibrated with analyses of national monitoring data (including non-linear relationships and ecologically meaningful breakpoints) and expert judgement. We used a novel expected index of desirability to quantify the model outputs. Secondly, we applied the BBN to three case study catchments in Ireland to demonstrate the implications of changes in stressor levels for ecosystem services in different settings. Four out of the seven significant relationships in data analyses were non-linear, highlighting that nonlinearity is common in ecosystems, but rarely considered in environmental modelling. Deficiency of riparian shading was identified as a prevalent and strong influence, which should be addressed to improve a broad range of societal benefits, particularly in the catchments where riparian shading is scarce. Sediment load had a lower influence on river Science of the Total Environment 851 (2022) 158065

Introduction
While not always explicitly stated, ecosystem services protection and restoration are central to many policies at the global (Sustainable Development Goals, Convention on Biological Diversity), European (Water Framework Directive, EU Green Deal and EU Biodiversity Strategy) and local scales (river basin management planning). Yet aquatic ecosystem services remain poorly integrated into water resources management decisions, including in freshwaters (AQUATAP-ES, 2020;Geneletti et al., 2020). This partly relates to the paucity of practical guidance on how best to link pressures and associated stressors to impacts on ecosystem services, a key step in ecosystem services assessment frameworks .
Rivers are a key part of the global hydrological cycle and are vital for the delivery of aquatic ecosystems services. They are hotspots of biological activity, supporting disproportionally high biodiversity compared to terrestrial ecosystems (Dudgeon et al., 2006), and producing a wide range of aquatic ecosystem services that benefit human societies for their survival, prosperity and well-being, such as the provision of water for abstraction and recreational uses, angling, and wildlife value (Feeley et al., 2017;Böck et al., 2018;Christie et al., 2021). The morphology and hydrology of river networks are widely manipulated by humans, which together with water pollution, mainly from sediment and excess nutrients, make rivers one of the world's most degraded and threatened ecosystems (Reid et al., 2019). Climate change is already adding to these pressures through changes in flow and temperature regimes (Gudmundsson et al., 2021). Stressors can interact (Davis et al., 2018;Birk et al., 2020;Lemm et al., 2021) and managing the complexity of their behaviour presents a major challenge (Hering et al., 2015).
Maintaining and improving the health of the aquatic environments in Europe is driven by the Water Framework Directive (EEA, 2018). Despite the considerable monitoring effort and management actions undertaken so far, only 40 % of Europe's surface waters are achieving satisfactory conditions (EEA, 2018). Ireland is no exception, with almost half of its rivers failing to meet their environmental objectives (EPA, 2021). Moreover, the number of high-status rivers in Ireland has declined precipitously from over 500 to 22 since the 1990s, and there have been recent signs of reversals in previous positive trends (EPA, 2020). This necessitates a paradigm shift in the approach to environmental management. Current river assessment focuses on a favourable ecological state based on biological components, supporting physico-chemical and hydro-morphological conditions, and on chemical pollutants (Poikane et al., 2014). These endpoints have been chosen because their conditions are of inherent interest and/or have established causal relationships with catchment pressures and associated stressors, but they are not necessarily well-suited to conveying the ecosystem services value of rivers. However, management measures must take account of this value because it places emphasis on objectives to which stakeholders, managers and politicians can relate (MEA, 2005;Kelly-Quinn et al., 2021).
Linking pressures with ecosystem services relies on a series of steps interconnecting abiotic and biological components, with an additional component of human value perception. These elements tend to be captured and represented using a range of different methods. Abiotic processes are frequently described quantitatively by physical-chemical (process-based) models, ecological dependencies tend to be restricted to empirical (datadriven) models or expert knowledge (either qualitative or quantitative), whereas socio-economic information is frequently only presented qualitatively or quantitatively using monetary values. Bringing such diverse information sources together in a common cause-effect chain is challenging (AQUATAP-ES, 2021). Bayesian belief networks (BBN) can handle this complexity (Smith et al., 2018). BBNs model relationships using conditional probabilities, i.e., the probabilities of each of the different outcomes of a state or process, under the different conditions of its influencing variables. The causal dependencies of each process can be then linked in a feedforward network to form a relationship chain between the input and output variables of interest using Bayes' theorem (Kjaerulff and Madsen, 2013;Feld et al., 2020). Importantly, BBNs can combine qualitative and quantitative information sources to inform the conditional probabilities in the same network. Because of this versatility, BBNs have been increasingly used in ecosystem services frameworks (Haines-Young, 2011;McVittie et al., 2015;Mantyka-Pringle et al., 2017;Forio et al., 2020).
Here, our overarching aim was to develop a BBN model for integration of ecosystem services considerations into river management that would ultimately form the 'backbone' of a decision support tool (the ProgRES tool; Kelly-Quinn et al., in press). To achieve that, we advanced an existing proof-of-principle BBN  to causally link catchment stressors with ecosystem services by (i) strengthening the evidence using weighted evidence from an expert workshop, and (ii) calibrating the model using national monitoring data and expert judgement. Relative to the starting BBN, we extended the range of ecosystem services in the BBN for a more integrated assessment (Grizzetti et al., 2016) and used individual environmental conditions as inputs instead of explicit management options to permit a greater flexibility. We then applied the improved BBN to three case study catchments in Ireland to demonstrate the implications of changes in the intensity of individual and combined stressors for ecosystem services and their trade-offs in different settings.

Material and methods
BBNs represent relationships among variables as a network of nodes interconnected with arrows, The direction, form and strength of each relationship are quantified through a set of probabilities. These give the conditional probabilities of occurrence of each possible categorical state of a target variable (child node) depending on the states of the variables (parent nodes) that directly influence it (Kjaerulff and Madsen, 2013;Scutari and Denis, 2014). Netica 6.09 was used to develop the BBN (Norsys Software, Vancouver, Canada), and GeNIe Modeler 3.0 was used for the strength of influence and sensitivity analyses (BayesFusion, Pittsburgh, PA, USA). GeNIe sensitivity analysis algorithm calculates a complete set of derivatives of the posterior probability distributions over the target node(s) over each of the numerical parameters of the network. These derivatives indicate the importance of precision of numerical parameters in the network for calculating the posterior probabilities of the target nodes (BayesFusion, 2020).

BBN inputs and outputs
The BBN underwent multiple iterations of development from the initial BBN . This was done, firstly, to encapsulate a broader selection of ecosystem services. Secondly, owing to concerns of the consulted experts about the advisability of using just a limited number of explicit land management options as the input variables in the BBN, these were replaced with explicit user control of each of the individual stressors that are predominant in the study area and can be managed realistically (Potschin and Haines-Young, 2016), and general river characteristics were used instead as inputs. This was because using a fixed range of management options as the inputs would limit the future applicability of the model, for instance if new measures were developed or when the effectiveness of a measure depended on the river characteristics, i.e. other exogenous variables. Thus, by having physico-chemical and other characteristics of the river as the direct input variables, a more general and flexible model could be produced. The intermediate nodes and connections were restructured to facilitate these changes.
After these iterations, the BBN inputs consisted of six environmental conditions which can be altered by catchment management: riparian shading level, sediment load, organic matter concentration (quantified by biochemical oxygen demand; BOD), phosphate concentration, oxidised nitrogen concentration (the sum of nitrate and nitrite) and total ammonia concentration (the sum of ionised and gaseous ammonia). The BBN inputs also included three general characteristics of the river (alkalinity, flow regime and coarse fish presence, the latter including species such as roach, pike and perch (Rutilus rutilus, Esox lucius and Perca fluviatilis, respectively), which are generally not manageable, but support applicability of the BBN in different catchment settings. And finally, the BBN included a climate scenario input (Table 1). Downscaled predictions of future mean temperatures are available for Ireland (Nolan and Flanagan, 2020). However, downscaled predictions of other climate change impacts, such as temperature extremes and flow regime impacts were insufficiently conclusive for Ireland to support robust implementation in the BBN (Fealy et al., 2018;Nolan and Flanagan, 2020). Thus, the current BBN likely underestimates the effect of future climate change.
The BBN has six output variables reflecting different provisioning, regulating and cultural ecosystem services or their biological/functional indicators, selected as the most likely to be valued by people and also the most likely to respond to catchment management. The number of mayfly species (insects of the taxonomic order Ephemeroptera) and the density of breeding pairs of white-throated dipper (Cinclus cinclus) were included as indicators of the value of wildlife as a cultural ecosystem service. Both positively correlate with river quality, and they appeal to a variety of river interest groups (Ormerod and Tyler, 1993;Ormerod et al., 1985;Kelly-Quinn and Bracken, 2000). Nutrient assimilation potential was included as a regulating service to encapsulate the capacity of the river to intercept phosphorus and nitrogen pollution, to self-purify and protect areas further downstream, including estuaries and marine waters. Recreational water quality (cultural service) and abstraction water quality (provisioning service) were included to reflect the potential for in-situ and ex-situ water uses, respectively. And finally, salmonid angling potential was included because it is a measure of a key recreational activity in Ireland that reflects a cultural ecosystem service (Table 2).

Expert knowledge
A one-day workshop was held on 29 January 2020 at University College Dublin with the participation of 17 experts from the national agencies with aquatic environment remits (e.g. Environmental Protection Agency, Inland Fisheries Ireland), and national and international universities (Table A.1 in Appendix A). The experts were asked to express their judgements on the nature and strength of the relationships between key nodes in the network by estimating the probability of responses conditional on selected environmental conditions (conditional probabilities). The experts were split into six teams of 2-3 participants, with each team providing independent judgements about the same set of network relationships that could not be determined independently from available data. In particular, each expert group was asked to estimate conditional probabilities for the following intermediate biological variables: 'Coarse fish density', 'Trout density', 'Salmon < 10 t km −2 y −1 Dod. Organic matter concentration (biochemical oxygen demand, mean; BOD) High density', 'Algal biofilm abundance' and 'Filamentous algae cover'. Probabilities were populated for each combination of the states of all parent nodes, so they captured any perceived interactions among them. Each group was asked to state their confidence level in their own expertise knowledge of each variable on a scale from 1 to 10. Probabilities across the six expert groups were then averaged after weighting by their stated confidence level for the BBN (Appendix A). The experts were also asked to comment on the BBN structure; in particular, if any intermediate variables or links were missing or redundant, and if the states of the BBN variables were adequately defined and quantified.

Other information sources
In Irish environmental law, Statutory Instrument (S.I.) number 77 of 2019 implements the European Union Environmental Objectives (Surface Waters) (Amendment) Regulations 2019 and sets boundary concentrations defining moderate, good and high status for water protection in Ireland, which were used here as the primary information source for defining model node states for oxygen, BOD, phosphate and ammonia. For these nodes, the legislative limits superseded any other information sources to link explicitly the BBN state specifications to the legal thresholds and facilitate their understanding and use by stakeholders. Because oxidised nitrogen does not have legislative thresholds, surrogate thresholds are used, based on the recommendations of the Environmental Protection Agency (EPA, 2011). Peer reviewed literature and reports from Ireland provided an evidence base for the biological effects of sediment in particular, either as a sole or co-occurring stressor (Garcia-Molinos and Donohue, 2009; Conroy et al., 2016;Davis et al., 2018). International sources were used throughout the BBN, but in particular to investigate the influences on benthic algae (Welch et al., 1988;Wagner et al., 2015), fish (e.g., Solomon and Lightfoot, 2008; with references therein), dipper (Ormerod and Tyler, 1993;Ormerod et al., 1985) and nutrient assimilation (Hall et al., 2009;Tank et al., 2018;Kikuchi et al., 2020). The findings of the SILTFLUX project  were used to inform sediment load thresholds.

National monitoring data
Data from the national river monitoring programme were a key source of quantitative information for calibrating the BBN. The data consisted of (i) counts of benthic macroinvertebrates from 2007 to 2018, collected once every three years from 176 sites; (ii) hydro-morphology (River Habitat Assessment Technique; RHAT) 2007-2018, usually collected together with the macroinvertebrate data; (iii) phytobenthos 2016-2018, collected once every year from 200 sites and (iv) physico-chemistry 2007-2018; collected at varying frequency from 4154 sites, but usually monthly from the biological sites.
The analyses comprised the examination of (i) data distribution, using histograms and summary statistics, (ii) linear and non-linear relationships between variables, using generalized additive models (GAM; Wood, 2017) and (iii) thresholds in individual data distributions or in the relationships between variables, using classification and regression tree analysis (CART; Hothorn et al., 2006;Appendix Table B.1). For the analysis of relationships and thresholds, physico-chemical data were screened for those locations matching the biological monitoring sites and aggregated (as the 5th percentile for dissolved oxygen and as the mean value for all else) over the time intervals of the biological data to match their temporal frequency.
GAM was used to initially analyse (i) how BOD, average summer water temperature and filamentous algae cover influenced dissolved oxygen saturation; (ii) how inorganic nutrient excess (quantified as the mean of total ammonia, oxidised nitrogen and phosphate, each first scaled by dividing by the Water Framework Directive pass/fail threshold), summer water temperature, shading score and siltation score influenced filamentous algae cover; (iii) how alkalinity, dissolved oxygen saturation and RHAT hydro-morphology score influenced mayfly richness; (iv) how shading score influenced summer water temperatures. Fewer data points were available for the variables filamentous algae and siltation score than for the other variables in the same model and thus these curtailed the model sample size. They were subsequently dropped in the final models of dissolved oxygen saturation and filamentous algae cover, respectively.
Each variable was Box-Cox transformed to reduce heteroscedasticity, then centred at its mean and scaled by its standard deviation. Repeated assessment rounds were included in each model as a random factor. GAM models were initially fitted with all terms as smoothers and replaced with linear terms one by one in the order of increasing estimated degrees of freedom. Non-linearity was detected for a particular explanatory variable, and a smoother term was retained, if a change of a smoother into a linear term increased Akaike information criterion by ≥2. In the other cases a linear term was chosen instead. We used CART on the untransformed data to seek break points in the same relationships. All data analyses were performed in R v4.0.2 (R Core Team, 2020).
These analyses were used to calibrate (i) numerical ranges of the different categories for a given variable (i.e. states), (ii) the relative strength of the different influences on a particular target variable (hierarchy), and (iii) the relative strength of the different intensities of a particular influence on a particular target variable (to identify non-linear relationships; Appendix B).

Expert judgement
The entire BBN was additionally 'sense-checked' by the project team using their expertise and scientific literature. The BBN was also checked to ensure that any adjustments of one parent variable had not distorted the effect of another, and if any indirect effect on a particular parent node did not distort the overall effect. This was done by examining the responses of posterior probabilities of child nodes to changes in the states of parent nodes. Furthermore, we checked if the overall behaviour of the network, including strength of influence and sensitivity, aligned with the general understanding of river functioning. In addition, the outcomes of the network applied to the three case study catchments were checked for alignment with the general knowledge of each catchment's functioning. These calibration steps prompted further rounds of adjustments. Finally, the BBN was tested using an online demonstration interface (Kelly-Quinn et al., in press) by external experts against their knowledge of the relevant cause-effect relationships, which resulted in one more adjustment refining the effects of increasing shading intensity.

Case study catchments
To demonstrate the BBN, we performed model simulations based on the catchment-scale data of selected three case study rivers representing  (2022) 158065 4 different characteristics and pressures ( Fig. 1): (i) River Dodder, a predominantly urban catchment in the east of Ireland; (ii) River Suir, an intensively agricultural catchment in the south of Ireland; and (iii) River Moy (part), a mixed land-use catchment with peaty soils in the west of Ireland. Owing to the size and high spatial heterogeneity of pressures in the Moy, we focused on one of its sub-catchments 'Moy SC 010' (area of 178 km 2 ), which has abundant coniferous forestry and rural housing with septic tanks. See Kelly-Quinn et al. (2020) for further description of these catchments. The BBN input states for these three catchments have been determined and set to 100 % probability based on 2018-2020 EPA data for water chemistry, aerial photography for riparian shading intensity (by the ESManage project), modelled values for sediment load (from equations developed by the SILTFLUX project; Bruen et al., 2017), Inland Fisheries Ireland data for coarse fish presence and Office of Public Works water level data for flow regime.

Generating BBN outputs
The average strength of the influences on each node throughout the entire BBN (based on Euclidean distance) and the sensitivity of the ecosystem service outputs to each of the environmental conditions inputs were analysed using GeNIe Modeler to investigate the general behaviour of the BBN, independently from any catchment settings (see BayesFusion, 2020 for details of these methods). For the three case study rivers, the effects of changing individual stressor intensities on individual ecosystem services were quantified based on the change in the probabilities of their states relative to the baseline conditions for each river. To investigate combined effects of stressors, we looked at binary interactions by descriptively comparing pairs of stressors. Thus, we compared the predicted improvements of ecosystem services from baseline conditions for each river resulting from simultaneous improvement of stressor pairs to their most desirable conditions with their combined individual effects. Greater simultaneous effects would denote synergism, i.e. stronger ecosystem-service benefit from simultaneous improvement of the two stressors than predicted from the sum of their individual effects, whereas smaller simultaneous effects would denote antagonism.

Expected index of desirability
BBN outcomes are typically presented as the change in the probability of a selected target state(s) of the output nodes in response to a change of inputs (e.g. Forio et al., 2015;Xue et al., 2017;Feld et al., 2020;Pham et al., 2021). When an output node has only two states, this approach captures all the available information. This is because probabilities across all states of a node always add to 1 (and their changes add to 0), so only n -1 states are free to vary. However, if an output node has more than two states, relative changes among multiple target or multiple non-target states cannot be completely determined from a change in the probability of an individual state, and yet they may be of interest. For example, in our BBN, the output node 'Mayfly richness' has four states: 'low', 'medium', 'high' and 'very high'. The last two are intuitively desirable, but if their cumulative probability change were selected as the target response, we would not be able to distinguish any changes from 'high' to 'very high', or from 'low' to 'medium', and yet they both would be ecologically meaningful. To reduce such information loss from the aggregation of changes to multiple states of a given output as the model result (but not information loss from discretisation also associated with BBN models), we quantified the posterior probabilities of each output node as an expected index of desirability. This was derived as weighted probabilities across all its states, with the weighting factor distributed at equal intervals between 0 and 1 from the least to the most desirable states, respectively (0, 0.5 and 1 for three-state nodes; 0, 0.33, 0.67 and 1 for four-state nodes). Our expected index of desirability is influenced by the probabilities of the output node states, and any changes to them, in a metric that intuitively ranges from 0 to 1 for any number of numerical or categorical states (corresponding to the full probability of the least and the most desirable states, respectively; Appendix C). Even though this index is constructed based on the probabilistic output of the model, it does not reflect a probability of an outcome, but rather the most probable outcome.

Results
The final BBN consisted of 36 nodes representing model variables, 76 links representing relationships between variables and 5018 conditional probability values quantifying these relationships ( Fig. 2; Appendices D and E).

Linear and non-linear relationships
Summer dissolved O 2 concentration 5 %-ile was found to decrease linearly with summer temperature (t = −3.93, P < 0.0001), and non-linearly with BOD (F = 67.33, P < 0.0001), whereby its decrease rate accelerated with increasing BOD (adjusted R 2 = 10.3 %). CART analysis identified BOD splitting values of 0.96 mg O 2 L −1 , 1.98 mg O 2 L −1 and 2.46 mg O 2 L −1 (Appendix B). Filamentous algae cover increased linearly with summer temperature (t = 5.88, P < 0.0001) and with nutrient concentration (t = 2.68, P = 0.0076; adjusted R 2 = 9.7 %). CART identified temperature splitting values of 13.8°C and 16.2°C. Mayfly richness changed nonlinearly with habitat quality (F = 3.88, P = 0.0161) and with alkalinity (F = 4.63, P = 0.0065), whereby it initially increased with both of these characteristics, but with no further increase beyond intermediate habitat quality, and with an apparent decline at high alkalinity (adjusted R 2 = 8.3 %). CART identified habitat quality RHAT score splitting value of 0.42. Summer water temperature decreased non-linearly with shading intensity (F = 10.64, P < 0.0001), whereby its decrease rate accelerated with increasing shading (adjusted R 2 = 24.8 %). CART identified shading score splitting value of 2, denoting light to moderate shading (Appendix B).

Sensitivity analysis of the general BBN model
Sensitivity analysis of the BBN showed that among the different input environmental conditions represented in the BBN, organic matter, closely followed by sediment load and then riparian shading deficiency had the strongest effects on the chosen ecosystem services (Fig. 3a). Organic matter pollution acts in the BBN in several ways, its decomposition lowers dissolved oxygen and contributes to excess nutrients, and it shares sources with Escherichia coli contamination (sewage and slurry; Fig. 2). Sediment, Fig. 2. Final BBN structure. Nodes in bold font are environmental condition inputs (left) and ecosystem service outputs (right). Arrow thickness denotes the strength of influence. Arrow colour denotes the sign of the correlation: green is positive, purple is negative, yellow is changing direction over the parent node gradient (unimodal). when deposited, smothers riverbed habitats and, when suspended in the flow, directly impacts water quality for abstraction and recreation. Riparian shading subdues filamentous algae growth and cools stream water, while promoting habitat quality and heterogeneity via associated inputs of foliage, coarse woody material, underwater roots and overhangs. Among inorganic nutrients, phosphate had the strongest effect on riverine ecosystem services, followed by total ammonia and oxidised nitrogen. All three contribute in our BBN to the nutrient excess and in addition, ammonia is toxic to animals in its un-ionised form (as dissolved gas, as opposed to ions). Future climate change was predicted to negatively affect ecosystem services through increased summer water temperature (Fig. 2). Among the different ecosystem services represented in the BBN, recreational water quality, closely followed by abstraction water quality, were the most sensitive to the input environmental condition, predominantly to sediment and organic matter (Figs. 2 and 3a). Salmonid angling potential was most sensitive to organic matter, followed by riparian shading deficiency and ammonia. Nutrient assimilation potential was most sensitive to riparian shading. The two indicators of wildlife value, mayfly richness and dipper density, had a similar sensitivity, strongest to organic matter.

Catchment-specific outcomes
All three case study catchments had high alkalinity, which is typical of many parts of Ireland. The Dodder and the Moy were flashy, whereas the Suir was slow flowing, thus the latter was more prone to sediment settling on the riverbed (Table 1; Fig. 3b, c). The Dodder and the Moy did not have coarse fish, whereas the Suir did. The latter's competition with salmonids, and to some degree also the salmonids' higher susceptibility to water quality deterioration owing to its slower flow, diminished the relative importance of ammonia in the Suir (Table 1; Fig. 3b, c).
The three catchments varied in the set of active stressors and their severity to some degree, but with some overlaps (Table 1). Increasing riparian shading was predicted to considerably benefit all ecosystem services, more so in the Suir, because of its high phosphate levels and slow flow. The latter is conducive to more time for heat transfer and thus, without shading, would lead to higher temperatures, and together with phosphate favour proliferation of nuisance algae, all of which can be inhibited by shading. The least benefit from increased shading was projected for the Dodder catchment, which had already a higher shading intensity than the Suir or Moy (Table 3). Decreasing sediment load was predicted to considerably improve various aspects of water quality in the Suir and Moy. Decreasing organic matter concentration could bring considerable improvements to all ecosystem services in the Moy. Reducing phosphate inputs was predicted to considerably benefit most ecosystem services in the Dodder and especially in the Suir which has the higher baseline concentration of phosphate. A reduction of total ammonia from the baseline levels was projected to benefit salmonid angling potential in the Dodder, but not so much in the Suir where multiple other stressors affect salmonids even when ammonia concentrations are low (e.g., phosphate, sediment, nitrate, low levels of shading and coarse fish). Decreasing oxidised nitrogen concentration had little influence in any of the catchments because it was never higher than 'medium' and phosphate is thought to be more limiting of algal productivity than inorganic nitrogen in freshwaters (Table 3).

Combined effects of stressors
Simultaneous reduction in stressor pairs from among oxidised nitrogen, total ammonia and phosphate always resulted in marginally stronger improvements of ecosystem services compared to the sum of their individual effects, possibly denoting weak synergism. The same was observed for the combination of sediment with either oxidised nitrogen or total ammonia (Table 4). The combined effects of sediment with organic matter or phosphate were usually synergistic. On the other hand, the combination of shading with organic matter or phosphate was always antagonistic. The simultaneous effects of shading combined with either oxidised Fig. 3. Sensitivity of ecosystem services to environmental conditions in the BBN represented by the thickness of the connecting shapes in the overall BBN (a) and for the three case study catchments (b, c). The height of the boxes denotes the overall relative influence of each environmental condition, and the overall relative sensitivity of each ecosystem service (in the same order in all three panels), across all scenarios represented in the BBN. nitrogen, total ammonia or sediment relative to the sum of their individual effects varied among rivers and ecosystem services between additive, synergistic or antagonistic (Table 4).

BBN outcomes and management implications
This study set out to develop the model as a tool to advance efforts to incorporate ecosystem services into decision making related to the protection and management of freshwater resources. The anthropocentric focus of the ecosystem services approach may help to communicate the urgency of addressing water quality declines (MEA, 2005). By linking these changes to anthropogenic stressors, our BBN can help identify potential beneficial or detrimental effects of management measures targeting these stressors.
Based on the aggregated information and BBN modelling, deficiency in riparian shading was identified in this study as a prevalent and highly influential stressor, which should be addressed to improve a broad range of river ecosystem services and thus societal benefits. This is particularly important to safeguard river biodiversity and functioning and would also help mitigate the effects of heatwaves and further climatic changes (Seavy et al., 2009;O'Briain et al., 2017). Notably, well-managed riparian vegetation can support stream functioning in multiple ways (Feld et al., 2018) and can generate societal benefits of their own, such as recreational and wildlife Table 3 Condition of ecosystem services in three river catchments expressed as an expected index of desirability ranging from 0 to 1 (i.e., worst to best; in bold) and predicted changes to these services resulting from alterations of environmental conditions. The direction and degree of change are indicated by the colour and intensity of shading: purple denotes a deterioration, green denotes an improvement.

RIVER MOY (part) Index
Mayfly richness (wildlife)  Table 4 Predicted improvements of ecosystem services from baseline conditions in three river catchments (measured as an expected index of desirability ranging from 0 to 1) resulting from improvement of stressor pairs to their most desirable conditions expressed as the difference between joint and individual effect of pairs of active stressors. Positive values denote synergism (green shading; stronger ecosystem-service benefit from simultaneous improvement of the two stressors than predicted from the sum of their individual effects), negative values denote antagonism (purple shading; weaker ecosystem-service benefit from simultaneous improvement of the two stressors than predicted from the sum of their individual effects). values (Nóbrega et al., 2020;Riis et al., 2020). Whereas oxidised nitrogen had little influence in any of the case study catchments, it controls algal productivity in coastal waters (O'Boyle et al., 2015;Malone and Newton, 2020), where it can affect a broad range of marine-based ecosystem services (Kermagoret et al., 2019), which were outside of the scope of this study. Integrated catchment management needs to consider such farreaching influences. The prevalence of the individual land-use-related stressors (sediment, organic matter, inorganic nutrients), and the ecosystem service benefits of their reduction varied among the three case study catchments to some degree, depending on the river conditions and interactions between stressors. Increasing riparian vegetation would be most beneficial where it is currently scarce (the Suir and the Moy catchments) and its ecosystem service benefits are likely to be more substantial in slow-than in fastflowing rivers (the Suir). Sediment load interacted synergistically with organic matter and phosphate in the Moy and the Suir, the two key drivers of eutrophication in freshwaters, producing a greater impact than expected from the sum of their individual influences across all ecosystem services in our BBN. This is likely because each of these pairs has a different and multifaceted mode of action. They are among the most ubiquitous stressors on rivers (Birk et al., 2020) and tackling them simultaneously was predicted by the BBN to yield additional societal benefits. Sediment load was also involved in some synergistic interactions with shading intensity, but this varied among rivers and ecosystem services.
Stressor rank identified in this study represents average conditions across the entire range of scenarios represented in the BBN, but it may vary according to site-specific conditions. In particular, inorganic nutrients (oxidised nitrogen, total ammonia and phosphate) ranked among the least influential environmental condition across all scenarios because their effect was partly conflated with organic matter, which both contributes nutrients through mineralisation, and has an additional direct effect on oxygen concentration. Thus, inorganic nutrients may be of higher importance in the absence of organic pollution. Another example of synergy is that riparian shading can have an even stronger effect on ecosystem services under high nutrient concentrations and future climate warming, because of the increased role of shading for cooling the stream water and subduing light for algal productivity. Ammonia in its toxic form can result in pulses of fish kills, but it is less likely to be in its toxic form in low alkalinity waters (IPCS, 1986), sediment is likely to get periodically washed out in flashy rivers , and so on. Furthermore, our BBN represents average conditions over the modelled catchment, but this can vary among specific locations, and thus a smaller spatial scale of modelling may be more appropriate for some applications, depending on the scales being considered. Therefore, specific management advice should always explicitly consider such context-dependencies and spatial heterogeneity. These findings may also inform the hypotheses of research seeking to disentangle multiple stressor effects/interactions, including direct and indirect effects on ecosystem services.
Each element of our BBN model underwent multiple iteration to assure that it is robust, but data were only available for the subset of the relationships. We trust the model predictions qualitatively, in terms of stressor importance hierarchy and the sensitivity rank of ecosystem services. However, the quantitative predictions of the likelihood or degree of improvement need to be interpreted with caution.

BBN methodological advances
Unlike the approach of only quantifying BBN model output based on selected target state(s) that is common in the literature (Forio et al., 2015;Xue et al., 2017;Feld et al., 2020), our expected index of desirability reduces information loss regardless of the number of node states. Although we used equal intervals of the weighting factor among node states, the numerical node states themselves were distributed at varying intervals (for 'Mayfly richness' and 'Dipper density'), reflecting their non-linear scaling with environmental conditions. Furthermore, the weighting factor among node states can be set at varying intervals if this better suits the objectives of a particular modelling scenario. The full response of any node with numerical states can be also captured by its expected value, derived as a sum of the product of each state and its corresponding probability , but this was not possible for our nodes with categorical states. Thus, our expected index of desirability is a simple, intuitive and flexible approach to capture the full BBN outcomes, for both numerical and categorical node states.
The Water Framework Directive requires comprehensive monitoring of aquatic ecosystems across the EU (EEA, 2018), and these data can be used to supplement expert opinions and other information sources for the BBN, as we have done here. Importantly, four out of the seven significant relationships were non-linear (Appendix B). This highlights that non-linearity is common in real ecosystems, but it is not considered in many ecosystem service models. We hope that our work will promote more targeted research on non-linearity and synergies in future decision support tools. Where comprehensive data exist, they can strengthen the BBN, provided that they are compatible with the BBN structure. A BBN can in fact be built entirely from data (Scutari and Denis, 2014;Feld et al., 2020), but this is not always possible, due to data limitations. The purpose of our BBN was to specifically connect ecosystem service outcomes to specific stressors, which dictated the BBN structure. This is not necessarily aligned with the goal of national environmental monitoring, which focuses on readily quantifiable abiotic and biotic conditions pertaining to ecological health per se. Such quantifiable metrics for the specific ecosystem services are yet to be developed and could be included in the national monitoring programmes. For example, we used mayfly richness as an indicator of wildlife value. The Irish monitoring data for river invertebrates are generally collected at a coarse taxonomic resolution, but it is reasonably fine for mayflies . Thus, our mayfly richness index could be approximated from the existing monitoring data and, consequently, we used it in data analyses. However, the available monitoring data was not adequate to quantify to the same extent any of the other ecosystem services.

Limitations
Broad-scale field data tend to contain multiple sources of variance, which may not all be quantified and thus are difficult to account for in statistical analysis. Consequently, even for well-known relationships, our analysis resulted in a relatively low coefficient of determination. Indeed, field data tend to be 'noisy', which complicates predictive modelling. Nonetheless, data analyses informed the strength and shape of some of the relationships, strengthening our BBN.
Expert opinions tend to be a key part of BBN development (McVittie et al., 2015;Constantinou et al., 2016). The conditional probabilities inherently incorporate uncertainty of outcomes by distributing them among the different node states, but this does not necessarily account for other sources of uncertainty, such as variability among experts. In our workshop, we split the expert pool into six groups to address this problem. The additional benefit of such an approach is that it down-weights a potential influence from particularly dominant personalities, which then becomes restricted to their own sub-group. This element itself is not novel (Kuhnert et al., 2010;Kelly-Quinn et al., 2020), but our weighting of expert opinions advances this method. We asked the experts to declare the degree of their confidence in their expertise relating to each node examined at the workshop and used this information in weighting their probability estimates. However, people can be overly conservative or overly confident so this can be improved in future through more objective metric of expertise (e.g. a publication output metric) or experience (e.g. number of years active professionally), and by validating the model against empirical evidence/experiments.

Conclusions
Management needs to reconcile the competing interests of landscape elements set against the costs and benefits of management actions. An ecosystem services modelling framework such as our BBN-informed decision support system can support this.
Key advances of the BBN methodology presented here include refined expert knowledge, analyses of national monitoring datasets and the expected index of desirability to quantify the outcomes. Relative to our starting BBN, we extended the range of ecosystem services. We incorporated a regulating service (nutrient assimilation potential), albeit only qualitatively owing to a knowledge gap in this respect. The BBN can be further refined by wider testing and emerging evidence.
While the model was demonstrated for three case study catchments, it can be generalized to some degree. It can be readily parametrised for most other Irish freshwater catchments by choosing appropriate input settings and possibly in the temperate zones more broadly by adjusting its structure. Furthermore, a similar general approach can be adapted to other regions, to other ecosystem types and to address different societal values. Developing and advancing support systems for complex decisionmaking scenarios, as we have done here, is important to inform management and policy.

Data availability
The authors do not have permission to share data.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.