Water Quality Sustainability Evaluation under Uncertainty: A Multi-Scenario Analysis Based on Bayesian Networks

Abstract

With increasing evidence of climate change affecting the quality of water resources, there is the need to assess the potential impacts of future climate change scenarios on water systems to ensure their long-term sustainability. The study assesses the uncertainty in the hydrological responses of the Zero river basin (northern Italy) generated by the adoption of an ensemble of climate projections from 10 different combinations of a global climate model (GCM)–regional climate model (RCM) under two emission scenarios (representative concentration pathways (RCPs) 4.5 and 8.5). Bayesian networks (BNs) are used to analyze the projected changes in nutrient loadings (NO₃, NH₄, PO₄) in mid- (2041–2070) and long-term (2071–2100) periods with respect to the baseline (1983–2012). BN outputs show good confidence that, across considered scenarios and periods, nutrient loadings will increase, especially during autumn and winter seasons. Most models agree in projecting a high probability of an increase in nutrient loadings with respect to current conditions. In summer and spring, instead, the large variability between different GCM–RCM results makes it impossible to identify a univocal direction of change. Results suggest that adaptive water resource planning should be based on multi-model ensemble approaches as they are particularly useful for narrowing the spectrum of plausible impacts and uncertainties on water resources.

1. Introduction

The maintenance of good water quality resources is essential to protect both ecosystems and human health, and they represent one of the main targets of both the European Water Framework Directive (2000/60/CE) and Sustainable Development Goals (i.e., SDG6) [1].

Changes in the global climate system are expected to have major consequences on the qualitative aspect of available water resources [2,3,4,5,6,7]; thus, assessing the impacts of future climate change scenarios on water systems is necessary to ensure a sustainable management of water supply for multiple purposes.

The inherent complexity, variability, and randomness of water systems, their interaction with socio-economic factors including the land use and population growth, and the high degree of uncertainty stemming from climate change make the assessment of climate change impacts on water resources particularly challenging [8].

Uncertainty plays a prominent role in climate change science and climate change impact science, with hydrology and water resources research in particular [8,9,10,11]. According to Parker et al. [12] and Hawkins et al. [13], it can be attributed to a number of reasons including (i) scenario uncertainty, arising from our limited understanding about the path of greenhouse gasses emissions and socio-economic development; (ii) internal climate variability, due to the inherent variability of the climate system components, processes, and their interaction; (iii) model uncertainty, caused by the different formulations used to represent climatic processes in climate and impact models.

A proper understanding of the type, sources, and effects of uncertainty is needed to achieve the goals of reliability and sustainability in water system management and planning under changing conditions [14,15]. Uncertainty quantification is vital to facilitate a risk-based approach to decision-making, where the range of possible futures is considered [16,17], and costs and benefits of adaptation are estimated accordingly. For this reason, uncertainties should be communicated as an inevitable component of each climate impact assessment study in a form which is also understandable by a non-scientific community to avoid misjudged information and to prevent overconfidence in impact projections [18].

A promising way to evaluate and deal with uncertainty is represented by ensemble modeling approaches [19]. Multi-model ensembles are commonly used to investigate structural model uncertainty, employing more than one climate model to perform multiple simulations and analyzing how climate change projections differ. The development of ensembles in both climate and impact studies is strongly encouraged by the Intergovernmental Panel on Climate Change (IPCC) since the Fourth Assessment Report (AR4, 2007) [20], which suggests the use of multiple climate models and scenarios to cover different sources of uncertainty [20]. The variability among ensemble components, in fact, can be used as a measure of the state of the knowledge. Furthermore, it could be useful to describe the confidence about the impact of climate change on the system modeled, supporting more robust decisions. In other words, if most ensemble members give comparable results, high confidence in projected climate change impacts is obtained, while, by contrast, if a large spread between components exist, there is less confidence in the forecasted impacts. Furthermore, it was shown that ensembles often give a more accurate prediction of future climate impacts than even the best individual model [21,22,23].

Relying on the extensive experience acquired in climate modeling, the use of ensembles was also transferred to the water resources field where attempts to build ensembles of impact models and scenarios (i.e., hydrological, water quality) are becoming increasingly common [24,25,26] to support water system management and adaptation.

In this respect, the paper proposes a Bayesian network (BN)-based approach to develop an ensemble of impact scenarios simulating the effect of different climate change projections on the quality of water of the Zero river basin (Italy). Accordingly, BNs are used as a modeling framework to evaluate the uncertainty due to global climate model (GCM)–regional climate model (RCM) structure and representative concentration pathways (RCPs), helping in determining and communicating the level of confidence of projected water quality alterations between baseline and future climate regimes.

BN outcomes (i.e., multiple impact scenarios) can be used to inform the spectrum of plausible effects of expected climate change on the Zero river basin and, thus, support the choice of effective adaptation strategies for a sustainable management of water resource quality at the local scale.

After a brief introduction to the study area (Section 1), this paper describes the methodology and input data employed (Section 2) and, finally, discusses the scenarios developed for the Zero river basin case study (Section 3), together with their uncertainty analysis.

Study Area

The Zero river basin (ZRB) (latitudes 45°28’ north (N)–45°48 N, longitudes 11°54’ east (E)–12°25’ E) (Figure 1) is located within the Venetian floodplain (northern Italy) and it is a sub-basin of the Venice lagoon watershed (Figure 1a), covering an area of 140 km². The Zero river (Figure 1b), which is 47 km long, together with the Dese rivers, provides the greatest contribution of freshwater (21% of the total) to the lagoon of Venice [27].

The basin features a Mediterranean climate but typical traits of more continental climates [28], with an average annual precipitation of around 1000 mm (period 2007–2012) and an average annual temperature of 14 °C (period 2004–2013). It is characterized by marked inter-annual climate variability, which can originate years climatologically very different from each other. The land use of the ZRB is mainly characterized by agricultural areas, representing 73% of the total surface, while the remaining surface of the basin is covered by artificial (24%), semi-natural, and forested areas (4%).

Agricultural areas are dominated by industrial crop typologies, including corn (45%) (i.e., Zea mays L.), soy (9%) (i.e., Glycine max L.), and autumn–winter cereals (13%) such as winter wheat (i.e., Triticum aestivum L.) and barley (i.e., Hordeum vulgare L.). A negligible percentage of the agricultural land is also used for the cultivation of beets and other permanent horticultural crops.

Artificial surfaces are mainly represented by housing areas (54%), industrial businesses (32%), and transportation and services (14%). Accordingly, several industrial and residential activities exist on the basin. Three wastewater treatment plans (i.e., Morgano, Zero-Branco, and Castelfranco Veneto) (Figure 1b) with capacities ranging from 2500 to 32,000 of population equivalents (PE) directly discharge into the Zero river.

The intensive agriculture, characterized by an elevated level of fertilization, and the dense urbanization represent significant pollution sources for the area and the main factor responsible for the excessive nutrient loadings in the reaching bodies of the Venice lagoon. Nutrient pollution is a major concern in the area considering the risk of eutrophication and toxic algae blooms which can threaten the good qualitative status of waters with consequent implications for environmental and human health [29]. Climate change, inducing extreme changes in temperature and precipitation trends, could exacerbate such nutrient pollution, altering those hydrological processes (e.g., runoff, river flow, water retention time, evapotranspiration) that regulate the mobilization of nutrients from land to inland and coastal water bodies.

2. Materials and Methods

To understand the complexity of the interactions between climate, hydrology and water quality parameters in the Zero river basin, we adopted the BN-based integrated approach proposed in Figure 2.

The core of the proposed approach is represented by a Bayesian network, which is employed as a modeling tool for the simulation of multiple nutrient loadings scenarios, and for the analysis of their uncertainty. Different data coming from hydrological models, historical observations, and expert judgment (

Table A1. List of input data used for the implementation of the BN in the Zero river basin.

Data Type	Description	Time Scale	Resolution	Source
Observations
Land cover map	▪ Land-use map of the Veneto region	2006	1:10,000	Regione del Veneto Infrastruttura dati territoriali (http://idt.regione.veneto.it/app/metacatalog/)
Climatic data	▪ Daily precipitation ▪ Maximum/minimum daily temperature ▪ Daily evapotranspiration	2004–2013	3 stations (i.e., Castelfranco, Veneto, Zero-Branco, Mogliano Veneto) (Figure 1)	ARPAV Servizio Meteorologico
Water quantity and quality data	▪ Observed daily river discharge ▪ Observed nutrient (NO3−, NH4+, PO43−) concentrations in the lagoon	2007–2012	2 stations (i.e., manual station (Code 122), automatic station (Code: B2q) (Figure 1)	ARPAV Servizio Acque Interne MAV (Magistrato Acque Venezia)
Point-source pollution	▪ Monthly N and P loadings from WWTP and industrial discharges	2004–2013	3 stations (i.e., Morgano, Zero-Branco, Castelfranco Veneto) (Figure 1)	ARPAV Servizio Acque Interne
Hydrological simulations
Water quantity and quality data	▪ Simulated runoff ▪ Simulated N and P load in the runoff	2004–2013	River basin	SWAT (Soil Water Assessment Tool) simulations [30]
	▪ Simulated river discharge Simulated nutrient loadings (NO3−, NH4+, PO43−) in the lagoon	2004–2013	1 station (i.e., manual station (Code 122)

, Appendix) were structured and combined into a probabilistic form to develop and train the BN at different levels of implementation.

Qualitative information elicited from experts was used to develop the conceptual model of the network and to train socio-economic and agronomic variables (i.e., irrigation, fertilizer application) of the model for which quantitative data were not available. Observations regarding the main climatic parameters (i.e., precipitation, temperature, and evapotranspiration) and point-source pollution sources (i.e., wastewater treatment plants (WWTPs) and industrial discharges), together with Soil and Water Assessment Tool (SWAT) model simulations (i.e., runoff, river discharge (Q)) for the current conditions developed by Reference [30], were used for the training of the network for the period 2004–2013. Additional observed hydrologic data (i.e., river discharge (Q), nutrient concentrations (i.e., NO₃⁻, NH₄⁺, PO₄³⁻)) coming from the water quality monitoring station were used to evaluate the performance of the model under current conditions.

After the training and validation, an ensemble of climate change projections generated by coupling different combinations of Global Climate Models (GCMs) with regional climate models (RCMs) was used as input for scenarios analysis to assess the effect (and uncertainty) of future climate change on nutrient loadings.

2.1. Climate Change Projections

To assess the effect of climate change on nutrient loadings (i.e., NO₃⁻, NH₄⁺, PO₄³⁻), changes in temperature and precipitation over future scenarios were selected as climate change indicators and used as input for the development of alternative nutrient loading scenarios using the BN model. The main aim of the study, however, was to capture uncertainties across a range of available GCM–RCM structures and representative concentration pathways (RCPs); thus, in order to represent the widest range of temperature and precipitation changes projected for the case study area, different climate change model outputs were considered (

Table 1. Global climate model (GCM)–regional climate model (RCM) projections selected and implemented in the Bayesian network (BN). SMHI—Swedish Meteorological and Hydrological Institute; DMI—Danish Meteorological Institute; CMCC—Centro Euro-Mediterraneo sui Cambiamenti Climatici.

No.	Global Climate Model (GCM)	Regional Climate Model (RCM)	Representative Concentration Pathways (RCPs)	Resolution	Time Range	Institute
1	HadGEM2-ES	RCA4	4.5, 8.5	12 km	1970–2099	SMHI
2	IPSL-CM5A-MR	RCA5	4.5, 8.6	12 km	1970–2100	SMHI
3	CNRM-CM5	RCA6	4.5, 8.7	12 km	1970–2100	SMHI
4	EC-EARTH	RCA7	4.5, 8.8	12 km	1970–2100	SMHI
5	MPI-ESM-LR	RCA8	4.5, 8.9	12 km	1970–2100	SMHI
6	CNRM-CM5	CCLM	4.5, 8.10	12 km	1950–2100	CLMcom
7	CMCC-CM	COSMO-CLM	4.5, 8.11	8 km	1976–2100	CMCC
8	HadGEM2-ES	RACMO22E	4.5, 8.12	12 km	1950–2099	KNMI
9	EC-EARTH	HIRHAM5	4.5, 8.13	12 km	1951–2100	DMI
10	EC-EARTH	RACMO22E	4.5, 8.14	12 km	1950–2100	KNMI

). This allowed including both the “worst” and “best” case scenarios in the BN, giving the users substantial flexibility in exploring and understanding the possible implications of climate change in the future. GCM–RCM combinations were selected among those available considering different kinds of features, including (i) their representativeness for the case study area and for the selected time periods (i.e., 2041–2070 and 2071–2100); (ii) their ability to perform at high spatial resolution; (iii) their availability in an open-source format.

Ensembles of 10 climate change models were selected (Table 1), including the CMCC-CM/COSMO-CLM GCM–RCM and nine GCM–RCM model combinations from the EURO-CORDEX project [31].

The CMCC-CM global model [32] is a coupled atmosphere–ocean general circulation model, while the COSMO-CLM (CCLM) [33] is a high-resolution (between 1 and 50 km) climate regional model; both were developed by the Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC), and, when coupled, they allow a spatial resolution of 8 km for the selected region.

EURO-CORDEX is the European branch of the CORDEX initiative sponsored by the World Climate Research Program (WRCP) with the aim of organizing an internationally coordinated framework to produce improved regional climate change projections for all land regions worldwide based on dynamical statistical downscaling models forced by multiple GCMs. CORDEX results were used as input for climate change impact and adaptation studies within the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC). In this study, nine climate change scenarios resulting from different combinations of GCMs and RCMs at 12-km spatial resolution were selected (Table 1). Different GCMs and RCMs were developed by different research groups including the Danish Meteorological Institute (DMI), the Swedish Meteorological and Hydrological Institute (SMHI), and the Met Office Hadley Centre (MOHC) (Table 1).

Based on the outputs of the selected GCM–RCMs (Table 1), different climate change scenarios were developed for the Zero river basin case study by extrapolating the mean temperature (°C) and the cumulative precipitation (mm) calculated on a monthly basis. Specifically, for each GCM–RCM, five different 30-year scenarios were developed for a control period (i.e., 1983–2012), a mid-term period (i.e., 2041–2070) and a long-term period (i.e., 2071–2100) under two different representative concentration pathways (i.e., RCP4.5 and RCP8.5). RCP4.5 represents the moderate emission scenario which predicts an increase in radiative forcing up to 4.5 W∙m⁻² by 2100 and a stabilization of the emissions (i.e., 650 ppm) shortly after 2100 [34], while RCP8.5 was chosen as representative of the extreme emission scenario, in which the greenhouse gas (GHG) emissions and concentrations increase considerably over the 21st century, leading to a radiative forcing of 8.5 W∙m⁻² by 2100 [35], thus describing a future without any specific climate mitigation target.

To make the outputs of GCM–RCMs suitable to be implemented at the spatial scale of impact assessment models, a bias correction was applied [30]. GCMs, in fact, have a spatial resolution too coarse for local-scale assessments and, for this reason, they are generally coupled with RCMs to consider the effects of orography, land–sea surface contrast, and land surface characteristics. However, RCMs also often show significant biases due to imperfect conceptualization, discretization, and spatial averaging within grid cells; therefore, a bias correction is required.

For the data used in this study, the linear scaling (LS) method was applied to correct the biases in the monthly values of temperature and precipitation based on observed ones. The LS method was applied using the software CLIME, a geographic information system (GIS) software for climate data analysis developed by the Regional Models and Geo-Hydrogeological Impacts (REMHI) division of CMCC, as extensively described in Reference [30]. Specifically, the method was implemented for all 10 climate scenarios for every weather station of the case study (Figure 1) using the rainfall and temperature observations for the period 1993–2012 as a correction factor. Once corrected, outputs of the GCM–RCMs for each of the 10 climate scenarios and for each of the three weather stations of the case study (i.e. Castelfranco Veneto, Zero-Branco, Mogliano Veneto) were elaborated to obtain suitable inputs for the BN model.

2.2. Bayesian Network Model

A BN model was employed to assess and compare the impacts of different climate change scenarios on nutrient loadings (i.e., NO₃⁻, NH₄⁺, PO₄³⁻) in the transitional waters of the Zero river basin, thus generating an ensemble of impact scenarios supporting the identification of climate change effect on water quality.

The BN was implemented by building on a BN model previously developed and validated in a case study [36] which was extended to allow the incorporation of multiple GCM–RCM inputs. The BN for the Zero river basin was developed and run using the software HUGIN Expert, version 8 [37,38]. For additional details about the methodology and data used to develop the BN, please refer to Reference [36].

2.2.1. BN Development and Training

The BN structure was designed following the DPSIR (Drivers-Pressures-State-Impacts-Responses) framework, starting from the conceptual model described in Reference [36]. An influence (i.e., “box and arrow”) diagram was developed including the most relevant systems variables (i.e., nodes), as well as the links between them (i.e., directed arcs), allowing the identification of the main cause–effect pathways between input variables, represented by climatic changes and land use, and output variables, represented by the increase in nutrient loadings (i.e., NO₃⁻, NH₄⁺, PO₄³⁻) discharged by the Zero river basin into the Venice lagoon. Successively, the BN was trained, assigning states, prior information and conditional probabilities to all nodes of the network, translating the conceptual model into an operative probabilistic form.

The training was performed using a heterogeneous set of information for the period 2004–2013 (Table A1 and

Table A2. Overview of nodes and states in the Bayesian network model for the Zero river basin.

Node	Description	Type	States	Parametrization Method
Season	Alternative seasons	Labeled	Winter	Expert judgement
			Spring
			Summer
			Autumn
Climate change scenario	Alternative climate change scenarios	Labeled	Baseline 1983–2012	CMCC-CM/COSMO-CLM simulations
			RCP 4.5 2041–2070
			RCP 4.5 2071–2100
			RCP 8.5 2041–2070
			RCP 8.5 2071–2100
Agricultural land scenario	Extension of land (ha) occupied by agricultural activities under different scenarios	Labeled	Actual 2004–2013;	Observations, LUISA simulations
			Future 2050
Temperature	Seasonal average temperature (°C)	Numeric interval	0–8.37	Observations
			8.37–13.79
			13.79–19.21
			>19.21
Precipitation	Seasonal cumulative precipitation (mm)	Numeric interval	0–201.50;	Observations
			201.50–328.73
			328.73–455.96
			>455.96
Potential Evapotranspiration	Seasonal cumulative potential evapotranspiration (mm)	Numeric interval	0–133.85	Observations
			133.85–228.3
			228.3–322.75
			>322.75
Effective rainfall	Seasonal cumulative effective rainfall reaching the soil (mm)	Numeric interval	0–64.13	SWAT simulations
			64.13–122.95
			122.95–181.77
			>181.77
Crop water needs	Seasonal water demand for different crop typology (mm)	Numeric interval	0–109.77	Equation [36]
			109.77–213.64
			213.64–317.50
			>317.50
Irrigation	Seasonal amount of water applied as irrigation	Numeric interval	<−55.29	Equation [36]
			−55.29–101.28
			101.28–257.86
			>257.86
N fertilizer application	Nitrogen fertilizer applied for each season according to different crop typology (kg/ha)	Numeric interval	0–45.74	Expert judgment
			45.74–87.52
			87.52–129.30
			>129.30
P fertilizer application	Phosphorus fertilizer applied for each season according to different crop typology (kg/ha)	Numeric interval	0–25.41	Expert judgment
			25.41–50.83
			50.83–76.25
			>76.25
N diffuse sources	Seasonal amount of nitrogen coming from agricultural practices (kg)	Numeric interval	0–7388.86	Equation [36]
			7388.86–13,959.99
			13,959.99–20,531.11
			>20,531.11
P diffuse sources	Seasonal amount of phosphorus coming from agricultural practices (kg)	Numeric interval	0–5169.28	Equation [36]
			5169.28–10,221.75
			10,221.75–15,274.21
			>15,274.21
N point sources	Seasonal amount of nitrogen coming from point sources (i.e., wastewater treatment plants (WWTPs) and industrial discharges) (kg)	Numeric interval	0–9382.64	Observations
			9382.64–10,389.82
			10,389.82–11,396.99
			>11,396.99
P point sources	Seasonal amount of phosphorus coming from point sources (i.e., WWTPs and industrial discharges) (kg)	Numeric interval	0–1143.64	Observations
			1143.64–1478.99
			1478.99–1814.35
			>1814.35
River discharge	Seasonal average river discharge (L/s)	Numeric interval	0–1458.96	SWAT simulations
			1458.96–2360.53
			2360.535–3262.102
			>3262.10
Runoff	Seasonal cumulative runoff (mm)	Numeric interval	0–49.90	SWAT simulations
			49.90–90.15
			90.15–130.40
			>130.40
N in runoff	Seasonal amount of nitrogen loaded in the runoff (kg/ha)	Numeric interval	0–0.63	SWAT simulations
			0.63–1.19
			1.19–1.75
			>1.75
P in runoff	Seasonal amount of phosphorus loaded in the runoff (kg/ha)	Numeric interval	0–0.44	SWAT simulations
			0.44–0.87
			0.87–1.30
			>1.30
Total N loading	Seasonal nitrogen load in the river (kg)	Numeric interval	0–17,031.20	Equation [36]
			17,031.20–24,401.92
			24,401.92–31,772.64
			> 31,772.64
Total P loading	Seasonal phosphorus load in the river (kg)	Numeric interval	0–5405.76	Equation [36]
			5405.76–9710.91
			9710.91–14,016.07
			>14,016.07
Loading NO3− lagoon	Seasonal loading of NO3− reaching the lagoon (kg)	Numeric interval	0–28,047.50	SWAT simulations
			28,047.50–48,615.00
			48,615.00–69,182.50
			>69,182.50
Loading NH4+ lagoon	Seasonal loading of NH4+ reaching the lagoon (kg)	Numeric interval	0–3224.52	SWAT simulations
			3224.52–5009.3
			5009.3–6794.17
			>6794.17
Loading PO43− lagoon (kg)	Seasonal loading of PO43− reaching the lagoon (kg)	Numeric interval	0–1978.90	SWAT simulations
			1978.90–2954.00
			2954.00–3929.10
			>3929.10

, Appendix) at seasonal time steps including historical observations, hydrological model simulations (i.e., SWAT), and expert opinions. Specifically, for nodes associated with climatic variables (i.e., temperature, precipitation, evapotranspiration), probabilities were learned directly from the frequency of observations of weather monitoring stations available in the case study (Figure 1). Nodes related with point pollution sources were trained using the nutrient loadings measured in the outflow from three different WWTPs (i.e., Morgano, Zero-Branco, Castelfranco Veneto) (Figure 1) located in the basin, summing up their respective contribution.

Probability distributions of hydrological variables (i.e., runoff, river flow, nutrients loadings, N and P in the runoff) were instead calculated based on the frequencies of results of hydrological simulations performed with the Soil and Water Assessment Tool (SWAT) [39]. Finally, nodes describing agronomic practices (i.e., water needs, irrigation, P and N fertilizer application) were trained through expert elicitation or by applying empirical equations due to the lack of quantitative information and experiences in the case study. An exhaustive description of assumptions and information used to train the BN can be found in Reference [36]. Figure 3 shows the configuration of the BN for the Zero river basin once states, prior information, and conditional probabilities of each node were parametrized.

2.2.2. BN Evaluation

The developed BN was evaluated with the main aim to assess if it was able to correctly represent nutrient loadings and dynamics in the case study. Specifically, two main forms of model evaluation were performed, namely, model predictive accuracy and sensitivity analysis [36].

Predictive Accuracy

Predictive accuracy assessment was performed by comparing BN simulations with observations from water quality monitoring stations available for the case study. Specifically, observed nutrient loadings were calculated, multiplying the observed water flow (Q) and nutrient concentrations measured at two different water quality monitoring stations located in the case study (Figure 1) managed by ARPAV Servizio Acque Interne and the former MAV (Magistrato alle Acque di Venezia). The manual station 122 (45°33’ N and longitude 12°15’ E), provided seasonal data, and the automatic station B2q (45°34’ N and longitude 12°17’ E) provided daily data. Both stations were specifically identified and they are routinely used for (i) the assessment of the good environmental status of the Zero river according to the requirements of the Water Framework Directive (WFD); (ii) the assessment of the compliance with the maximum admissible load of nutrients discharged into the Venice lagoon from the drainage basin fixed by the national competent law (DM 09/02/1999).

For this reason, the two stations were considered particularly representative to measure the condition of the Zero river and its basin, thus providing reliable data for the evaluation of BN performance.

Specifically, loadings of nitrogen nitrate (NO₃⁻) and ammonium (NH₄⁺) were obtained using data from station B2q, while loadings of phosphate (PO₄³⁻) were calculated using data from station 122. A consistent set of observations was available only from year 2007 to 2012 and, therefore, the evaluation was conducted only for this period.

For each output node (i.e., NO₃⁻, NH₄⁺, and PO₄³⁻ loadings), correctly classified instances (CCIs) were assessed as the percentage of cases correctly predicted divided by the total number of cases, providing the measure of how many instances the model predicted correctly when tested against known case outcomes (i.e., observations). Error rates, used as evaluation criteria, were then computed and depicted in confusion matrices as suggested in Reference [40].

Sensitivity Analysis

Another form of evaluating the developed model entails sensitivity analysis, which allows testing the sensitivity of model outcomes to variations of model parameters [41]. Sensitivity analysis in BNs can measure the sensitivity of outcome probabilities to changes in input nodes or other model parameters, such as changes in node type of state and their coarseness; therefore, it useful to detect the most relevant variables within the network. Sensitivity analysis was performed using two types of measures: entropy and Shannon’s measure of mutual information [42].

The entropy (H(x)) of the probability distribution of a variable (x) expresses the measure of the associated uncertainty of the random process with a particular probability distribution (P(x)) [42]; it is calculated using the following function:

H(x) = −∑ P(x) log P(x).

(1)

Reducing entropy by collecting information, in addition to the current knowledge about the variable x, is interpreted as reducing the uncertainty about the true state of x. Accordingly, the entropy function enables an assessment of the level of uncertainty/certainty about the state of the output node and of the additional information required to specify a particular alternative.

Entropy can be seen as a score of a variable richness (i.e., how much information is within the data for that particular variable) [43,44], and it was used to rank nodes from the most uncertain to the least uncertain, where the most uncertain variables are the least informative within the network.

In addition, the sensitivity of one node to multiple other nodes was evaluated using Shannon’s measure of mutual information (MI) as follows:

MI (Y,X) = H(Y) − H(Y|X).

(2)

MI enables assessing the effect of collecting information about one variable (Y) on reducing the total uncertainty about variable X. When MI is equal to zero, the condition of one node does not affect the state of the other and, therefore, the nodes can be defined as mutually independent [43].

2.2.3. Scenarios and Uncertainty Analysis

The model developed as above was used in this study to perform scenario analysis, allowing the assessment of the relative change in outcome probabilities of nutrients under different climate change conditions (Section 2.1), thus obtaining an ensemble of multiple impact scenarios. For each GCM–RCM combination (Table 1) and climate change scenario (i.e., 2041–2070 and 2071–2100 under two different representative concentration pathways, RCP4.5 and RCP8.5), the probability distributions of temperature and precipitation were calculated based on the frequency in the respective model simulations (Section 2.1). The BN was then run, alternatively fixing the evidence of being in a particular scenario by assigning 100% probability to the related state in the “climate change scenario node”, letting the information propagate through nodes linked by conditional probability (Figure 3) and calculating the change in the posteriori probabilities of output variables (i.e., NO₃⁻, NH₄⁺, PO₄³⁻ loadings).

Moreover, uncertainties in projected loadings due to the application of the ensemble of 10 GCM–RCM couples was performed by comparing outputs obtained with each of the different GCM–RCM combinations across scenarios and seasons. Specifically, the changes in the probability of each loading class in the mid-term (i.e., 2041–2070) and long-term (i.e., 2071–2100) simulated periods were compared against the respective baseline scenario (i.e., 1983–2012) for each combination of GCM–RCM.

3. Results

3.1. BN Evaluation

3.1.1. Accuracy

As described above, a data-based evaluation was performed to assess the ability of the model to correctly predict instances in an independent dataset. Accordingly, BN predictions were tested against observations from water quality monitoring stations available from ARPAV for the case study (Figure 1), generating confusion matrices representing the percentage of CCIs and, consequently, the error rates (Figure 4). Observations were available only for 2007–2012 and, therefore, the evaluation was conducted only for this period.

In addition, the expected values of the probability distributions of nutrient loadings (i.e., NO₃⁻, NH₄⁺, PO₄³⁻) for observed data were compared with those obtained through the Bayesian network outputs (Figure 5).

Overall, the BN was able to reproduce the observed nutrient dynamics with loadings closely replicated for most seasons. The evaluation produced very good results for phosphate (PO₄³⁻), while, for ammonium (NH₄⁺) and nitrate (NO₃⁻), the correlation between observed and predicted nutrient loadings was slightly worse. Indeed, overall, the BN was able to correctly classify 87.50% of instances for PO₄³⁻, 63.64% for NH₄⁺, and the 66.67% for NO₃⁻, when tested against the observed dataset (Figure 4). The BN overpredicted the decrease in ammonium and nitrate loadings between spring and summer, while it underestimated the autumn loadings (Figure 5) for all three nutrient species (i.e., PO₄³⁻, NH₄⁺, NO₃⁻) and the winter loadings of NH₄⁺ and NO₃⁻.

3.1.2. Sensitivity Analysis

Entropy was calculated for each node of the network (

Table 2. Node ranking according to entropy score.

Variable	Entropy H(x)
Effective rainfall	3.26
Season	1.39
Temperature	1.38
N point sources	1.36
River flow	1.34
N runoff	1.33
N diffuse sources	1.33
Total N Loading	1.31
Evapotranspiration	1.28
Loading NO3−	1.24
Irrigation	1.24
P point sources	1.17
Runoff	1.12
Precipitation	1.06
N fertilizer application	1.04
Loading NH4+	0.98
Loading PO43−	0.97
P runoff	0.96
P diffuse sources	0.96
Water needs	0.56
P fertilizer application	0.56

), allowing us to rank all variables from the most uncertain to the least uncertain, where the most uncertain variables were those characterized by high entropy and were, thus, the least informative probability distributions.

Results showed that the variable characterized by the least informative probability distribution and a particularly high value of entropy was “effective rainfall” (3.26); however, variables characterized by intermediate values of uncertainty were nodes directly influencing “loading NO₃⁻”, including “N point sources” (1.36), “river flow” (1.34), “N runoff”(1.33), “N diffuse sources” (1.31), and “total N Loading” (1.31). These uncertainties propagated through the network and, as a result, “loading NO₃⁻” was the output node characterized by the highest value of entropy (1.24), while, for others, the uncertainty was moderate (0.98 and 0.97).

Table 3. Summary of the mutual information (MI) analysis presenting the top five most influential variables on output nodes.

Sensitive Node	Node Affecting Sensitivity	MI
Loading NO3−	River flow	0.54
	Total N loading	0.30
	N diffuse sources	0.23
	N runoff	0.23
	Evapotranspiration	0.23
Loading NH4+	River flow	0.37
	Loading NO3−	0.18
	Total N loading	0.17
	Runoff	0.15
	N diffuse sources	0.11
Loading PO43−	Total P loading	0.53
	P runoff	0.39
	P diffuse sources	0.39
	River flow	0.34
	Runoff	0.31

provides a ranking of the top five most influential variables on output nodes based on the mutual information analysis. The output nodes “NO₃⁻ loadings” and “NH₄⁺ loadings” were both highly sensitive to “river flow” (MI = 0.54 and 0.37, respectively). “PO₄³⁻ loadings” resulted highly sensitive to “total P Loading” (M = 0.53), “P runoff” (M = 0.39), and “P diffuse sources” (M = 0.39).

In general, hydrological variables, which in turn were strongly influenced by climatic ones (i.e., precipitation, temperature), were those most influential on other network variables. By contrast, variables related to agronomic practices and land use had a mild effect (MI < 0.1) on other variables with the exception of “N fertilizer application”, which moderately affected “NO₃⁻ loadings” (MI = 0.22). In particular, point sources had negligible effects on all output nodes (MI < 0.04) with respect to diffuse sources.

3.2. Climate Change Scenarios for the Zero River Basin

Different climate change scenarios were developed for the Zero river basin case study by extrapolating the mean temperature (°C) and the cumulative precipitation (mm) calculated on a monthly basis, based on the outputs of the selected GCM–RCMs (Table 1). Specifically, for each GCM–RCM, five different 30-year scenarios were developed for a control period (i.e., 1983–2012), a mid-term period (i.e., 2041–2070), and a long-term period (i.e., 2071–2100) under two different representative concentration pathways (i.e., RCP4.5 and RCP8.5).

Figure 6 shows the variability of temperature for different time periods and RCPs across different climate change scenarios used to inform the BN. It is possible to observe that the temperature variability across the future projection was quite narrow.

All climate scenarios agreed on projected temperature during the control period (i.e., 1983–2012). Greater variability, instead, was depicted for RCP8.5 where one model in particular (i.e., MPI-ESM-LR/RCA4, Model 5) of the ensemble projected lower temperatures in spring and higher temperatures in autumn. In general, all models predicted an increase in mean seasonal temperature with respect to baseline across the different considered scenarios (Figure A3, Appendix). MPI-ESM-LR/RCA4 (Model 5) represented the only exception, predicting a decrease in temperature in spring for RCP8.5 (Figure A3, Appendix). A greater increase in temperature with respect to baseline was predicted by RCP8.5 for the period 2071–2100.

By contrast, precipitation featured marked variability in all scenarios, as shown in Figure 7. All 10 GCM/RCMs of the ensemble generated quite similar statistics for the control period (i.e., 1983–2012) with a narrow range between maximum and minimum values for all seasons. By contrast, the variability increased consistently along the century, especially for RCP8.5. Greater variability can be seen in summer, autumn, and winter, where the range between maximum and minimum values projected by different GCM/RCMs became quite wide. However, while for winter and autumn most models agreed on an increase in the cumulative precipitation (Figure A4, Appendix), for spring and summer, models gave the opposite results, making it impossible to agree on the direction of change (i.e., decrease or increase).

3.3. Hydrological Responses to Climate Change

The BN model was run, alternatively fixing the probability distribution of precipitation and temperature according to the medium- and long-term projections (i.e., 2041–2070, 2071–2100) provided by the different available combinations of GCM–RCMs under two different representative concentration pathways (i.e., RCP4.5 and RCP8.5). Accordingly, the network was used to develop multiple impact scenarios linking the effect of future climate change projections on nutrient loadings. The developed scenarios represent the probability of different classes of nutrient loadings (i.e., low, medium, high, very high) calculated by the BN model as a result of changes in the probability distribution of input variables (i.e., temperature and precipitation).

Figure 8 gives a concise overview of the probabilistic results obtained through the BN for each season and scenario across the different GCM–RCM models considered (Table 1). Specifically, each triangular portion of the graph represents one of the different climate change scenarios considered (i.e., RCP4.5 2041–2070, RCP8.5 2041–2070, RCP4.5 2071–2100, RCP8.5 2071–2100), while, inside them, each slice represents the results of different GCM–RCMs arranged in a clockwise direction (i.e., from 1–10, in Table 1). Each slice, in turn, is divided into the four different classes of loadings with an amplitude corresponding to the value of the associate probability (i.e., from 0–100).

With regard to NO₃⁻ (Figure 8a), the impact scenarios reported that higher loadings will take place in autumn and winter, while the lowest loadings are predicted for summer. Across different models, in fact, in autumn, higher probabilities were associated with high (i.e., 48,615-69,182 kg/season, orange) and very high loading classes (i.e., >69,182 kg/season, red). The highest loading was predicted by the MPI-ESM-LR/RCA4 (Model 5) under the RCP8.5 2071–2100 scenario with 70% probability associated with the high loading class (Figure A5, Appendix).

In summer, by contrast, a higher probability was associated with low (i.e., 0–28,047 kg/season, green) loading classes, with the CMCC-CM/COSMO-CLM (Model 7) predicting the highest probability (77%) under the long-term RCP8.5 scenario (Figure A5, Appendix).

For ammonium (i.e., NH₄⁺), results across different models predicted high probabilities of low loading during summer and spring (Figure 8b). The lowest loading was predicted by the CMCC-CM/COSMO-CLM (Model 7) under RCP8.5 2041–2070 with a 97% probability of the low loading class (i.e., 0–3224 kg/season, green) (Figure A6, Appendix). In autumn, the probability of low loading states decreased gradually across the scenarios, followed by an increase in the probability of medium (i.e., 3224–5009 kg/season, yellow) and very high loadings (i.e., >6794 kg/season, red), respectively reaching 38% and 24% under RCP8.5 2071–2100 in the simulation with IPSL-CM5A-MR/RCA4 (Model 2).

Results for PO₄³⁻ showed a marked seasonality with high autumn loads and low loads in spring and summer across different scenarios (Figure 8c). In summer, in fact, higher probabilities were associated with the low loading state (i.e., 0–1978 kg/season, green). Specifically, the lowest loadings were predicted by the CMCC-CM/COSMO-CLM (Model 7) under the medium- and long-term RCP8.5 scenarios with a probability of 98% (Figure A5, Appendix). High loadings were instead predicted for autumn with probabilities of high (i.e., 2954–3929 kg/season, orange) and very high classes (i.e., >3929 kg/season, red) increasing across scenarios. The IPSL-CM5A-MR/RCA4 (Model 2), which described the most extreme loadings for the season, predicted probabilities of 34% and 16% of being in very high and high classes under the long-term RCP8.5 scenario (Figure A5, Appendix).

Uncertainty Analysis

The variability of results was also analyzed by comparing outputs obtained with each of the 10 GCM–RCM combinations across scenarios and seasons. To make the results comparable, the change in the probability of each loading class with respect to the baseline scenario (i.e., 1983–2012) was calculated for each combination of GCM–RCM.

Accordingly, in Figure 9, which provides an example for PO₄³⁻ loadings, negative values describe a decrease in probability of specific loading classes (i.e., colored bars) with respect to baseline, while positive values indicate an increase. Orthophosphate (i.e., PO₄³⁻) loadings showed clear variability during spring and summer. During these periods, in fact, half of considered models predicted an increase in loading, while others predicted a strong decrease. Less marked variability, however, was depicted under RCP8.5 2071–2100, where most models agreed on a reduction in loadings in the summer–spring period and an increase in probabilities associated with the low class. A good agreement among models, instead, can be depicted in autumn especially under RCP8.5 2071–2100, where most models predicted an increase in probabilities of very high and high loadings. Despite the good agreement on the increase in loading, a moderate variability in the magnitude of the change with respect to baseline remained. For RCP8.5 2041–2070, for instance, the maximum variation was related to MPI-ESM-LR/RCA4 (Model 5) (i.e., +20%), while EC-EARTH/RCA4 (Model 4) predicted an increase of +1.5%. In RCP8.5 2071–2100, the increase in probability ranged from +28% for IPSL-CM5A-MR/RCA4 (Model 2) to 2% for EC-EARTH/RCA4 and EC-EARTH/HIRHAM5 (Models 4 and 9). Also, in winter, a general increase in loading was predicted with an increase in probabilities associated with higher classes and a consequent decrease in probabilities of lower classes. The maximum increase (i.e., +10%) was depicted with EC-EARTH/RACMO22E (Model 10) under RCP4.5 2071–2100.

Results for NO₃⁻ and NH₄⁺ loadings presented a similar tendency (Figure A6, Appendix). The best agreement among models resulted for the autumn season, where an increase in loading was predicted across all scenarios and for all GCM–RCM combinations. Specifically, for NO₃⁻, an increase in probability of the high loading class was depicted, while, for NH₄⁺, the increase was associated with the highest loading class (i.e., very high). Also, for winter, the variability of results was quite low, with most models agreeing on an increase in probability of high and very high classes across different scenarios. By contrast, two models (i.e., CNRM-CM5/CCLM and CMCC-CM/COSMO-CLM (Model 6 and 7)) predicted a decrease in loadings for both NO₃⁻ and NH₄⁺. Large variability resulted for both summer and spring seasons; hence, it was not possible to identify a clear direction of change.

Overall, the results for different nutrient species highlighted that, in general, the best agreement between models resulted for autumn and winter, especially for RCP8.5 scenarios. In summer and spring, instead, variability was high and, thus, there was less confidence in the changes projected. This seasonal pattern of variability to some extent reflects that of precipitation (Section 3.1, Figure 6 and Figure 7), suggesting that this variable could play a major role in the model in determining both the direction and the magnitude of changes in nutrient loadings. Comparing the results (Figure 9, Figure A3 and Figure A4, Appendix) with the changes in precipitation across the different models (Figure 9 and Table A1 and Table A2, Appendix), a strong correlation between the increase in precipitation and increase in the probability of high loading can be found.

In summer and spring, in fact, those models which predicted the highest increase in probability of high loadings were also those showing a positive variation (i.e., increase) in precipitation with respect to baseline.

In this context, it is interesting to notice how, in spring (Figure A6 and Figure A7, Appendix) some models (e.g., Models 8–10) contemporarily predicted an increase in probability of the two most extreme classes (i.e., low and very high). The same models were also those predicting the highest increase in precipitation with respect to baseline (Figure A7, Appendix), suggesting that the unexpected high probabilities of very high loading classes could be related to the projections of extreme precipitation events in the considered scenarios.

4. Discussion and Conclusions

A BN was used to develop an ensemble of impact scenarios to investigate the potential mid- and long-term impacts of climate change on the nutrient loadings in the waters of the Zero river basin in northern Italy, one of the main tributaries of the Venice lagoon. Moreover, the uncertainty related to the implementation of 10 GCM–RCM combinations forced with RCP4.5 and RCP8.5 emission scenarios was analyzed.

The BN used was implemented by building on a model previously developed and tested in a case study [36], integrating a heterogeneous set of data coming from multiple information sources (i.e., observations, hydrological model simulations, climate change projections). The BN was evaluated through a cross comparison between predicted and observed loadings, providing satisfactory results on a seasonal time step; therefore, it is considered suitable for projecting future climate change scenarios.

According to the analysis of future climate for the Zero river basin, all projections agreed on an increase in the mean seasonal temperature with respect to baseline for both RCPs.

By contrast, precipitation featured marked variability across scenarios; while for winter and autumn most models agreed on an increase in the cumulative precipitation, for spring and summer, some models gave opposite results, making it impossible to agree on the direction of change (i.e., decrease or increase). The variability increased consistently along the century, especially for RCP8.5.

The impact scenarios developed showed that seasonal changes in precipitation and temperature are likely to affect nutrient loadings and, consequently, the water quality of the Zero river. Results suggest with good confidence that, across the considered scenarios, nutrient loadings will increase, especially during the autumn and winter seasons. Most models, in fact, agreed in projecting a high probability of an increase in nutrient loadings with respect to the current conditions. In summer and spring, instead, the large variability between different GCM–RCM results made it impossible to identify a clear direction of change.

The results were consistent with those obtained by Reference [45] applying the SWAT model for simulating the effect of climate change on hydrological and ecological parameters in the same case study. However, while conclusions for autumn are similar to those reached by other authors [26,46,47] for similar catchments in Europe and the United States (US), for spring and summer, the results differ. Xu et al., in particular, found that, in the Lake Erie region, spring loading of P will increase under RCP8.5 scenarios driven by an anticipation of snow-melting processes. Such discrepancies can be attributed to local and regional climatic characteristics which should, therefore, be taken into account carefully, as they can have a significant role in governing nutrient transport dynamics.

In the Zero river basin, nutrient loadings were found to be particularly sensitive to hydrological variables (i.e., river flow, runoff, N and P in runoff) directly correlated with climate variables (i.e., precipitation, temperature) and diffuse pollution, especially considering that most dramatic changes (e.g., increases in precipitation and runoff) will happen during seasons characterized by intensive agricultural activities (e.g., manure application, irrigation).

In spring and summer, in fact, NO₃⁻ and NH₄⁺ are commonly applied as fertilizers. In dry and warm conditions, NH₄⁺ is readily adsorbed to clay mineral and is, therefore, scarcely prone to movements; however, it becomes easily available in autumn, driven by runoff and extreme precipitation events. NO₃⁻, on the other hand, is highly soluble and, thus, suitable to be transported by hydrological flow. In autumn, the elevated temperature and wet conditions projected will enhance the nitrification process, making NO₃⁻ highly available. This, combined with the seasonal increase in the river flow, could explain the great increase in NO₃⁻ loading during autumn. In the soil, the soluble form of phosphorus (PO₄³⁻) is mobile, and it can be transported by diffusion or by surface water flow. At elevated temperatures and in dry conditions, however, PO₄³⁻ is easily adsorbed to clay particles or immobilized by organic matter accumulating in the upper soil layer. In autumn, an increase in runoff, following the enrichment of the topsoil of phosphorus occurring during the summer, increases PO₄³⁻ transport and, thus, its loading in the river. In addition, the projected increase of dry prolonged conditions in summer might speed up soil erosion phenomena and, consequently, enhance the runoff of adsorbed mineral forms of phosphorus through the basin, leading to peaks of PO₄³⁻ loading in autumn as soon as the drought breaks.

At the same time, the large uncertainty in spring and summer loadings makes it difficult to predict the possible implications for the trophic state of the Venice lagoon. A significant increase in spring and summer nutrients delivered by the Zero river, during the season of growth for most phytoplankton species, would significantly increase the risk of harmful algae blooms and eutrophication phenomena.

The BN was revealed to be a suitable tool to characterize and communicate uncertainty on the effect of climate change and land use on water quality attributes in a policy-relevant manner; however, it is important to also acknowledge some limitations. Some uncertainty exists mainly due to the availability and quality of input data, especially regarding agronomic practices. Due to data constraints, in fact, fertilizer application and irrigation were considered uniform across the whole catchment, while they could vary considerably, both spatially and temporally. Furthermore, due to scarce information regarding point pollution sources, nutrient (N and P) loading was considered while neglecting to take into account the type of WWTPs and how they work in cases with a large amount of inflow water, for instance, during extreme precipitation events.

Improving the accuracy of input data throughout the catchment and involving a higher number of experts in the model development would improve its calibration, validation, and results.

Finally, changes in land use (i.e., agricultural land extension, crop typology distribution) and agricultural management practices (i.e., amount of fertilizer application), which were kept constant over future scenarios in this BN version, should be accounted for in future model improvements to provide a realistic picture of future risks threatening water quality sustainability. Accordingly, further improvements of the proposed approach will consider the implementation of a dynamic version of the BN [48] to better handle temporal dynamics over future scenarios, while also integrating land-use change projections.

Overall, the results obtained from this study show that the selection of climate change information to feed impact studies should be considered carefully as it strongly affects the outcome and the conclusions of the assessment. Studies based on only one GCM–RCM combination should be interpreted with caution, as results are highly dependent on the assumptions of the selected combination. Adaptation and management decisions are taken based on this information with the consequence that societies may underprepare for real risks affecting water systems, increasing the likelihood of severe impacts, or, by contrast, they may overreact, wasting resources and efforts targeting irrelevant threats.

Accordingly, an adaptive water resource planning method should be based on ensembles and multi-model probabilistic approaches rather than on an individual scenario and a single-value projection for the future. Through the identification of worst- or best-case scenarios, it is possible to bound the spectrum of plausible climate change impacts into an uncertainty space, inside which a set of optimal adaptation strategies can be defined and tested for the sustainable and climate-proof management of the water system.