Ecosystem Services Evaluation of Nature-Based Solutions with the Help of Citizen Scientists

Abstract

Ecosystem services are increasingly being considered in decision-making with respect to mitigating future climate impacts. In this respect, there is a clear need to identify how nature-based solutions (NBS) can benefit specific ecosystem services, in particular within the complex spatial and temporal dynamics that characterize most river catchments. To capture these changes, ecosystem models require spatially explicit data that are often difficult to obtain for model development and validation. Citizen science allows for the participation of trained citizen volunteers in research or regulatory activities, resulting in increased data collection and increased participation of the general public in resource management. Despite the increasing experience in citizen science, these approaches have seldom been used in the modeling of provisioning ecosystem services. In the present study, we examined the temporal and spatial drivers in nutrient delivery in a major Italian river catchment and under different NBS scenarios. Information on climate, land use, soil and river conditions, as well as future climate scenarios, were used to explore future (2050) benefits of NBS on local and catchment scale nutrient loads and nutrient export. We estimate the benefits of a reduction in nitrogen and phosphorus export to the river and the receiving waters (Adriatic Sea) with respect to the costs associated with individual and combined NBS approaches related to river restoration and catchment reforestation.

1. Introduction

Ecosystem services (ESs), identified in relation to the benefits provided to people from ecosystems, are increasingly being considered in decision making by modeling specific services and their geographical context and extent [1,2,3,4]. As in any model, the outputs are sensitive to the characteristics of the input data [5]. Often the availability of sufficient data to both develop the model and validate it is a major challenge.

Citizen science (CS) allows for the participation of trained citizen volunteers in research or regulatory activities, often resulting in an increased data collection on spatial and temporal scales [6,7,8,9]. CS also facilitates increased participation of the general public in resource management [10,11,12,13]. Despite the diffuse experience in ecological and environmental projects [14,15,16], CS approaches have seldom been used in ES studies focused on provisioning services [17,18], compared to regulating and cultural services [19,20,21]. Likewise, data from CS have the potential to support reporting SDGs indicators [11,22,23], in particular those related to SDG 6, clean water and sanitation.

ESs delivered by freshwater ecosystems are strongly influenced by the interaction of land management, catchment hydrology and climate conditions. How this interaction plays out over the complex geological, hydrological and socio-economic conditions of a catchment has direct consequences on ESs related to water quantity and quality [24,25]. Nutrient retention is a key ES, which was modeled using a range of approaches [5,26] that balance the complexity of model accuracy, feasibility, computational cost and applicability. Models that are based on well-tested, relatively simple algorithms using readily available data are more easily understood and adopted by stakeholders [27] to develop more integrated catchment management, fundamental for both river functioning as well as that of its receiving waters. Coastal zones receiving river inputs are experiencing major eutrophication with related environmental and economic damage [28].

Extreme eutrophication in Europe [29,30,31], North America [32,33] and China [34,35] have had major consequences on local populations and economies. This is particularly important in Europe, where states have invested significantly to achieve good ecological status in water bodies, required by the Water Framework Directive [36] and the Common Agricultural Policy (CAP) [37].

The success of these policies, particularly in Europe [38,39,40], the United States [41,42] and China [43,44,45], led to a significant reduction in P load in many rivers and marine waters [46,47,48]. Because N loading, largely associated with diffuse agricultural sources, was not simultaneously reduced, a significant excess of N over P in river loading became the rule, and eutrophication problems persisted in coastal zones. Measures taken to reduce N losses from agriculture were less effective than those devoted to urban wastewater treatment, both due to the difficulties in tackling diffuse versus point sources of pollution as well as due to agricultural intensification in some areas [46]. However, recent studies integrating agricultural activities for the quantification of diffuse nutrient loading (N especially, and P) and nutrient delivery to the marine coastal zone showed success [49,50,51].

This led regional and national decision-makers to explore effective and verifiable measures to reduce nutrient loading through NBSs. These efforts are most often focused on improved land use, with particular attention to the riparian interface between terrestrial and aquatic ecosystems, which plays an important role in nutrients and sediments dynamics [52,53,54,55,56]. The riparian forest provides protection against erosion, retention of pollutants, excessive nutrients runoff and water temperature [57,58,59,60], with benefits for the quality of human life and biodiversity [61,62,63,64]. In this context, Zhang et al. [65] and Yang et al. [66] emphasized the importance of improving our understanding of the benefits of NBS focused on riparian zone restoration.

Besides agricultural intensification, climate change can also enhance eutrophication and its effects on the aquatic environment. Changes in precipitation will modify nutrient delivery through changing surface and river water flows [67,68]. These effects are sensitive to precipitation intensity as well as aggregate precipitation as reduced river water flow increases residence time and decreases nutrient dilution [69,70,71]. Ecosystems’ responses to both climate change and eutrophication may be hard to unravel and thus difficult to predict, but current knowledge suggests that climate change will worsen the present situation [28,72].

The present study utilizes multiple data sources and CS to explore the link between NBS and ESs related to nutrient retention. We compare different riparian zone NBS in relation to their ES costs and benefits related to nutrient retention under a future climate change scenario (2050) in one of Italy’s most important rivers.

2. Materials and Methods

2.1. Study Site

The Piave River, in the northeast of Italy, extends from the Eastern Alps to the Adriatic Sea. It is the fifth-largest Italian river with a watershed of 4127 km². Approximately 300,000 people live in the catchment, which extends across three Italian regions. The upper catchment (from 500 to 1800 m a.s.l.) is covered largely by evergreen and deciduous forests (mainly conifers and broad-leaved trees), while alpine pasture/prairies and rock emergence dominate the highest elevations (above 1800 m a.s.l.) and with limited agricultural activities (Figure 1). Agricultural land use is more prominent in the lower part of the catchment (below 500 m a.s.l) with an increase in impervious surfaces [73].

2.2. Model

Model development and validation were performed using information from different regional, national and global datasets and through in situ monitoring, both from the regional environmental agency (ARPAV) and citizen scientists. Geographic and temporal dynamics of the nutrient export and retention in the Piave River catchment were modeled following the mass balance approach of the InVEST nutrient delivery ratio model (version 3.9) [74]. The model is based on mapping nutrient sources and their potential for transport to the river to identify the spatial variation in nutrient retention across the watershed with respect to different land use/land cover conditions (LULC, e.g., vegetative areas) and catchment morphology.

Catchment nutrients dynamics were based on nutrient loads across the landscape and the nutrient retention capacity of the landscape. Nutrients loads per LULC were identified from data acquired in collaboration with the regional environmental agency, river consortiums and available empirical data. Nutrient flows were divided into sediment-bound (transported via surface flow) or dissolved (and transported via subsurface flow). A nutrient delivery ratio index (NDR) was simulated for each pixel (20 m resolution) of the catchment, based on the nutrient loads, LULC, and a digital elevation model. At the watershed outlet, the P and N export to the river was calculated based on the weighted aggregation of pixel-level contributions: site-specific information related to the maximum retention efficiency, a runoff proxy (i.e., annual precipitation) representing the spatial variability in runoff potential, and an estimate of the proportion of nutrients delivered via subsurface and surface flows. The subsurface flow was considered for nitrogen only and was estimated from nitrate concentrations in the ground and surface waters.

2.3. Data Sources for Model Development and Validation

Data required to develop and validate the nutrient N and P models were obtained from multiple sources and include:

• A digital elevation model (DEM), of 20 m resolution, obtained from the Italian ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale) [75], was corrected to fill hydrological sinks and checked with the digital watercourse network to ensure routing along the specific watercourse, using QGIS 3.12 (https://www.qgis.org, accessed on 24 March 2021);
• LULC raster data (2018) were obtained from Corine Land Use Land Cover IV Level for Italy [76], at a 100 m resolution;
• Nutrient runoff proxies were based on raster precipitation data from 2018, provided by the river catchment authorities (Autorità di Bacino Distrettuale delle Alpi Orientali) at a 20 m resolution. These raster data were interpolated using an inverse distance weighting of information from 72 stations;
• Future Nutrient runoff proxies were based on precipitation estimates for 2050 by Deszi e al. [77] with high-resolution gridded surfaces at 1 km cell size developed in an Albers Equal Area Conic projection for Europe. Vector delineation of the watershed was obtained from the geoportal of ARPAV, relative to the Water Protection Plan 2015 (Piano di Tutela delle Acque) [78];
• The threshold value for flow accumulation, the number of upriver cells that flow into a cell before it is considered part of a river, was set to 1000, after several tests to compare the river layer output of the model to the measured river network data [79].
• Borselli’s k for the connection of the surrounding land to the river with respect to the ratio of nutrients reaching the river was set to 2 [80];
• The nutrient (N and P) sources associated with each LULC class (kg ha⁻¹ y⁻¹) were based on 2001 data from ARPAV for Corine LULC classes 111–243 (artificial surfaces and agricultural areas) for each municipality and scaled for relative population changes in 2020. For Corine LULC classes 3112–523 (forest and semi-natural areas and water bodies), nutrient load data were obtained from the ARPAV relative to 2018 [81]. The proportion of subsurface N, a floating-point between 0 and 1, was obtained by intersecting a Corine LULC 2018 vector layer with a groundwater infiltration potential layer [82]. In order to obtain the proportion of subsurface N, the groundwater N infiltration risk potential was compared to the protective soil capacity layer [83,84] (

  Table 1. The proportion of subsurface N weighted for the % of Corine LULC per risk class.

  Soil Protective Capacity Table			Assigned Risk Potential	Proportion of Subsurface_N to Surface N per Risk Class
  Protection Capacity	Infiltration Flux	Loss of NO3	Risk	Reference Values (to Be Weighted per % of LULC Polygons)
  High	<12%	<5%	Very low	0.05
  Medium high	12–28%	5–10%	Low	0.075
  Medium low	29–40%	11–20%	Medium	0.155
  Low	>40%	>20%	High	0.2

  ). For each LULC, the final subsurface_N value is weighted for the % of Corine LULC polygons per risk class and calculated as the median of all level 1 Corine LULC classes (i.e., all urban classes, all agricultural classes, all forest classes, etc.) (Table 1).

The proportion of subsurface N value was estimated using groundwater nitrate data from ARPAV for 2018 [85], which has a lower temporal resolution (30%) with respect to surface NO₃.

• Retention efficiency for N and P, as the maximum nutrient retention expected from each LULC type, were calculated following Pärn et al. [86], Mayer et al. [87] and Zhang et al. [65];
• Retention lengths for N and P for each LULC class, as the typical distance necessary to reach the maximum retention efficiency, were based on previous studies of riparian buffers [65,87] and ranged from 10 to 300 m. In the absence of data, the retention length was set to the pixel size;
• Subsurface critical length, the distance after which the soil retains N at its maximum capacity, was set to 200 m following Mayer et al. [87]. Maximum retention of N reached through subsurface flow was set to 0.8 [87];
• Nitrate and phosphate concentrations were obtained through regulatory (ARPAV) and CS measurements. Quarterly ARPAV monitoring of total nitrogen and total phosphorus in four stations in the upper catchment and four sites in the lower catchment were used for model development. Data are available online at www.arpa.veneto.it/dati-ambientali/open-data/idrosfera/corsi-dacqua (accessed on 3 August 2021);
• Trained citizen scientists used the FreshWater Watch method (https://freshwaterwatch.thewaterhub.org/content/freshwater-watch-how-guide, accessed on 3 August 2021) to determine nitrate and phosphate in 12 sites in the lower catchment for model validation. Online and video training followed the standard training program of this global citizen science project [88]. Thirty-five participants in the lower catchment of Piave from Ponte di Piave to the sea collected both observation data (color, presence of algae, etc.) and semi-quantitative measurements of water quality (nitrate, phosphate, nephelometric turbidity) [89,90]. Nitrate (NO₃-N) and phosphate (PO₄-P) were measured in closed plastic tubes, which are designed to mix a fixed volume of water with reagents to produce increasing color values (peak absorption at 540 nm) with increasing concentration. PO₄-P concentrations were estimated colourimetrically using inosine enzymatic reactions in seven specific ranges from 0.02 mg L⁻¹ to 1.0 mg L⁻¹ PO₄-P (<0.02, 0.02–0.05, 0.05–0.1, 0.1–0.2, 0.2–0.5, 0.5–1.0, <1.0 mg/L) [91]. Nitrate–nitrogen concentrations were estimated colourimetrically using N-(1-napthyl)-ethylenediamine [92], in seven specific ranges from 0.2 mg L⁻¹ to 10 mg L⁻¹ NO₃-N (<0.2, 0.2–0.5, 0.5–1.0, 1.0–2.0, 2.0–5.0, 5–10, >10 mg/L). The median values for each classification were used to allow for quantitative analysis. During the measurements, geo-location and time were recorded automatically using the FreshWater Watch app and transferred to the online database after validation (https://freshwaterwatch.thewaterhub.org/, accessed on 3 august 2021). Once entered, all data underwent quality control by project leaders and citizen scientists.

The complete biophysical table is available in the Supplementary Materials (Table S1).

2.4. Climatic Variations

The potential of NBS to influence ESs is moderated by the local climate conditions, particularly the seasonal and interannual variations in precipitation dynamics. The potential influence of climate change on ES delivery in the present study was explored using historical time series and global climate models.

The climate of the study area has a temperate humid regime with months of maximum precipitation in autumn and late spring (May), with minimum precipitation generally recorded in February and July. The average annual rainfall is highly variable. An analysis of local trends in precipitation and temperature was performed using:

• ARPAV data of Ponte di Piave station (1995–2018) lower catchment, Santo Stefano di Cadore station (1998–2018) upper catchment [93];
• NOAA data (National Oceanographic and Atmospheric Centre, Silver Spring, Maryland, USA) from NCEP reanalysis (1948–2018) [94];
• Global climate projections (2041–2070) [77].

The climatic components most directly impacting ecosystem services related to the water, carbon and nutrient cycle are the average monthly precipitation, the number of days with precipitation, average monthly temperature and average monthly maximum temperature. The data were subjected to decomposition methods to explore seasonal, inter-annual trends and teleconnections with global climate indices. The North Atlantic Oscillation (NAO) is an important teleconnection scheme influencing European climates. It is based on a dipolar model of mean sea level pressure over the North Atlantic extending from subtropical to sub-arctic latitudes. It is associated with variations in westerly winds relative to Western Europe, an important factor for winter weather in Europe [95,96].

Future projections were developed with the delta method based on multi-model CMIP5 projections for the 2011–2040 and 2041–2070 periods (referred to as the 2020s and 2050s). The baseline reference data for monthly precipitation were constructed with parameter regression of independent slopes model (PRISM) by Daly et al. [97]. We selected a median emission scenario, represented by (RCP) 4.5 that globally predicts a +1.4 °C (±0.5) by 2050.

2.5. NBS Scenarios

Corine Land Cover (CLC) was used to conduct the analysis for the current scenario (2018), with a spatial resolution of 100 m [98]. For the 2050 scenario, three alternative NBS land cover scenarios were developed based on consultations with river authorities and an analysis of photometric and satellite-based images. This analysis showed a consistent and continuous increase in riparian vegetation within the Piave river corridor [99] over the last five decades and a more recent trend in agricultural abandonment in the upper catchment.

Alternative NBS scenarios were generated using QGIS (version 3.12). The first (B1) was based on the reforestation of selected agricultural areas of the upper catchment (conversion of CLC class 243 to class 313). The second (A1) was based on the creation of a riverine hygrophilous forest (class 3116) and an erodible area (class 331) along the river corridor of the lower catchment. The final scenario was the combination of both A1 and B1. River downscaling was performed to highlight better the areas that contribute most to the nutrient delivery to the river. Eight scenario combinations for 2018 and 2050 were explored.

2.6. Cost Analysis of N and P

In order to compare the costs and benefits of different NBS scenarios, an estimate of the economic value associated with phosphorus and nitrogen retention service of the Piave catchment was determined using the reported costs of recent NBS projects focused on nutrient reduction. The projects considered were performed in the same region by the Drainage Authority ‘Consorzio di Bonifica Acque Risorgive’, from 2003 to 2020 [100]. These projects were financed by the Decree of the Ministries of the Environment and of Labour to reduce the nutrient loads to the Venice Lagoon. We conducted a multiple linear regression after adjusting for inflation to estimate the typical costs of N and P removal. These were compared to recent studies on the costs for nitrate and phosphate removal by improvements to wastewater treatment [101,102]. In general, cost estimation of nutrient reduction is based on the analysis and estimation of different components such as energy consumption, chemicals’ consumption, personnel salaries, maintenance expenses, construction materials and their quantities, mechanical equipment and land cost [103,104]. As NBS actions include costs that are not directly related to individual components, using an empirical aggregated approach was deemed more appropriate.

3. Results

The changes in nutrient retention-related ESs in the Piave catchment showed the relative impact of different NBS scenarios and the impact of climate change, in particular precipitation.

3.1. Climate Change

In 2018, the highest annual precipitation (mm) occurred in the upper catchment with a mean of 1599 mm/year, while the precipitation mean in the lower catchment was 1070 mm/year. The wettest month was October (396 mm), the month with the least rain was December (10 mm) (Figure 2).

Teleconnections with major climate indices indicated potential links of the catchment climate to conditions in the rest of Europe. The NAO index shows an increased variability since 1950 (Figure 3), associated with the changes in circulation due to differences in pressure at sea level. There were several long periods in which the anomalous circulation of the NAO persisted, particularly from the 1940s to the 1970s, with a downward trend and lower than normal winter temperatures. After a period of relatively negative NAO values, the last 7 years (2012–2018) have seen a significant increase, associated with warmer-than-usual winter temperatures and changes in the precipitation regimes in much of Europe.

The monthly dynamics of precipitation showed a much lower seasonality with respect to temperature, with the former showing a strong autocorrelation from 1995 to 2018 (Figure 4 and Figure 5).

There were no significant interannual trends (α = 0.01) identified in the seasonal decomposition for temperature and precipitation. However, there was an increase in the irregular component for temperature (p < 0.001), indicating a potential increase of more than 1.5 °C by 2050. The increase in the monthly maximum is also significant (p < 0.001), characterized by an increase (1.2 °C by 2050). The irregular component of precipitation (mm and days for months) did not show a clear trend. A clear negative relationship between the monthly NAO indices and precipitation (monthly and precipitation days) was identified. However, the links between the NAO regimes and the hydroclimate are widely regarded as not being constant over time [105]. The output of climate models [106,107] suggests a moderate reduction in expected rainfall in the coming decades. The same models estimate an increase in the temperature of the study area, which implies greater evaporation and a greater hydrological deficit. Current and future trends and variability will have clear impacts on ecosystem services related to nutrient flow within the catchment.

Future projections represent an ensemble average of 15 Atmosphere–Ocean General Circulation Models of the CMIP5 multi-model data set, corresponding to the IPCC Assessment Report 5 [108]. We selected the median emission scenario RCP 4.5 that globally predicts a +1.4 °C (±0.5) by 2050 [77]. This scenario confirmed a decrease in precipitation in 2050, with an estimated minimum of 758 mm/year and a maximum of 1968 mm/year (Figure 6).

3.2. Nutrient Export

In order to validate the modeled N and P load exported to the river, the ARPAV and citizen scientist data were used, considering monthly averaged nutrient concentrations and river water flow for 2018 [109]. The monthly averaged flow rate at Ponte di Piave (

Table 2. Piave river flow recorded at Ponte di Piave by flow velocity probe to the year 2018.

Piave Flow m3/s	Min	Max	Mean	SD
Ponte di Piave	8.2	582.2	116	184.96

) was used, which did not differ significantly from measurements made at Nervesa della Battaglia (130 m³/s).

ARPAV samples (

Table 3. Elaboration of ARPAV data from Ponte di Piave to the year 2018.

ARPAV	Min	Max	Mean	SD	Export ARPAV [ton/year−1]
T-N mg/L	1.19	2.35	1.62	0.43	5.75 × 103
T-P mg/L	0.016	0.059	0.029	0.01	1.04 × 102

) show high concentrations of total nitrogen (TN) in winter of 2.17 mg/L and lower concentrations in summer and autumn of 1.32 mg/L. Total phosphorus concentrations (TP) were high in autumn (0.037 mg/L) and low in winter (0.021 mg/L).

Trained citizen scientists collected more than 100 samples along the river in the lower catchment (

Table 4. Elaboration of Citizen Science data to the year 2020.

CS	Min	Max	Mean	SD	Export CS [ton/year−1]
NO3-N mg/L	0.1	7.5	1.37	0.61	4.48 × 103
PO4-P mg/L	0.01	1.15	0.32	0.22	8.04 × 102

), where the number of environmental agency monitoring stations is limited. Nitrate–nitrogen (NO₃–N) concentrations were high in autumn (1.63 mg/L) and low in summer (0.90 mg/L). Phosphate–phosphorus (PO₄–P) concentrations were low in winter (0.11 mg/L) and high in summer (0.52 mg/L).

The NDR model outputs showed a total N and P export of nearly 4500 tons/year of nitrogen and 600 tons of phosphorus (

Table 5. Nutrient Delivery Ratio model results referred to the current scenario (2018).

Nutrient NDR Model	Export_Tot [ton/year−1]	Surface Load [ton/year−1]	Sub_Surface Load [ton/year−1]
T-N	4.47 × 103	2.99 × 104	5.93 × 103
T-P	5.96 × 102	1.14 × 104	-

), with surface loads dominating over subsurface loads.

By using the precipitation scenario raster for 2050, changes in nutrient export were explored for each scenario of land use, with, without and combined NBS. By applying the same land cover in 2050, P export is expected to increase by 1.7% (

Table 6. NDR model results comparing current (2018) and 2050 estimates of nutrient loads and exports (exp) for the land-use scenarios: A0 B0, without modification into the catchment; A0 B1, with reforestation on the upper catchment; A1 B0, with increased riparian vegetation and erodible bank area in the lower catchment; A1 B1, applying NBS approaches to both the upper and lower catchment.

Year	Surf_Load_P [ton/year−1]	Exp_Tot_P [ton/year−1]	Surf_Load_N [ton/year−1]	Sub_Load_N [ton/year−1]	Exp_Tot_N [ton/year−1]
Current	1.14 × 104	5.96 × 102	2.99 × 104	5.96 × 103	4.47 × 103
2050—A0, B0	1.13 × 104	6.07 × 102	2.91 × 104	5.83 × 103	4.39 × 103
2050—A0, B1	1.07 × 104	4.52 × 102	2.87 × 104	5.80 × 103	4.14 × 103
2050—A1, B0	1.12 × 104	6.04 × 102	2.91 × 104	5.82 × 103	4.39 × 103
2050—A1, B1	1.07 × 104	4.48 × 102	2.86 × 104	5.80 × 103	4.13 × 103

). A decrease in the delivery and retention of total N and total P reflected the precipitation reduction for the 2050 projection [77].

The reforestation scenario (B1) provided a reduction both for N and P load and export (Table 6). P export decreased by 24%, corresponding to a reduction of 144 tons/year. The N export reduction was relatively lower (7%), with a reduction of 328 tons/year. The reforestation scenario was based on the conversion of CLC class 243 (Agriculture with natural vegetation) to class 313 (Mixed Forest) (11,549 ha) and 3116 (Hygrophilous Forest) (3025 ha) (Figure 7a), resulting in an increase of 5% in forest cover in the upper catchment

For the lower catchment riverine corridor scenario (A1), there was an estimated increase in the P export by 1.3% due to the increment of sand and gravel where the P is more linked. N export decreased by 1.8%, or 80 tons/year (Figure 8). This scenario was based on an increase in the hygrophilous forest along the lower catchment river corridor from 6.54% to 34.87% (class 3116 area 650 ha) and the erodible area between the river banks (class 331 area 420 ha) (Figure 7b).

When both NBS scenarios were considered (A1, B1), there was a decrease in both export and load for both nutrients (Figure 8). N export decreased by nearly 8%, or 340 tons/year, while P export decreased by 25%, 148 tons/year (Figure 8).

4. Discussion

4.1. Nutrient Dynamics and Distribution

There is growing awareness that NBS can help to reduce climate change-related impacts as well as providing ecosystem services [110,111]. The results of the simulated nutrient conditions in four different NBS scenarios for 2050 suggest that the nutrient load and export in the Piave river catchment is sensitive to climate change and to an increase in wooded areas, particularly in agricultural areas where there is a trend in agricultural land abandonment [112,113]. This rural exodus, which began after World War II, triggered the process of natural vegetation regrowth that continues into the present [114]. Reforestation presents multiple benefits, promoting biodiversity [115] as well as reducing surface water runoff and soil erosion [116], control sediment loss and improve soil properties [117].

The overall impact of each NBS scenario is clear in the reduced export of both phosphorus and nitrogen, with the largest reduction for a single scenario occurring for B1, the reforested area of the upper catchment (export reduction for N of 328 tons/year, export reduction for P of 144 tons/year). The values confirm support studies that show increased forest cover and decreased agricultural areas decrease sediment, nitrogen and phosphorus exports [55].

The reductions associated with the scenario of A1 in the lower section are less evident but still represent a significant reduction in nutrient export. The lower catchment is heavily impacted by intensive agriculture both inside and outside the river banks, leaving very little area for hygrophilous riparian vegetation, with that remaining colonized often by invasive species. NBS actions focused on river renaturalization and restoration were shown to provide multiple benefits on nutrient mitigation as well as other ES. The erodible bank expansion in the A1 scenario is also intended to restore river dynamics but, as a consequence, will increase P export [118].

While the overall reductions in nutrient export provide a tool for understanding the benefits of NBS related reforestation and restoration activities, key information can be gained from exploring the spatial distribution of ecosystem services based on different scenarios. By using spatially explicit models, it was possible to identify the locations within the catchment with the greatest sensitivity to climate change and NBS actions, in this case with respect to nutrient loads and nutrient export. As expected, the areas closest to the river network have the highest nutrient export in relation to their retention capacity of the vegetation downslope, substrate and drainage.

The difference between current and 2050 A1B1 scenarios shows that areas closest to the river have the largest potential for a reduction in P export (Figure 9). Changes in N export (Figure 10) were less spatially distinct as nitrate is more mobile [110]. Complete 2050 scenario maps are available in the Supplementary Materials (Figures S1 and S2).

4.2. ESs Evaluation of NBS

The selection of the best single or combined NBS depends on the relative value associated with each ecosystem service and the increase or decrease in that service provider. By considering only nutrient reduction, it is possible to compare the impact of each scenario and considering overall changes in main drivers related to climate and land use. The results based on the simulated nutrient load and delivery under four different NBS scenarios, A0B0, A1B0, A0B1 and A1B1, show significant differences. There is a clear potential for a reduction in nitrogen export in 2050 for all scenarios. This was associated with reduced expected precipitation of 1363 mm/year compared to 1653 mm/year in 2018, particularly in the elevated alpine parts of the catchment. In fact, the largest reduction in nitrogen export was achieved by reforestation of the upper catchment (B1). This NBS alone provides a potential reduction of 328 tons/year. The costs of natural reforestation are limited, and the benefits to lower catchment and receiving waters (Adriatic Sea) could be translated to support for farmers to manage these lands as productive forests. It is expected that the new European agricultural policy will incentivize this transformation of land use. On the other hand, N and P loads from urban wastewater represent a significant source [119]. Recent studies on lowland rivers that have created riparian buffer strips have shown an elevated efficiency to remove nitrate [56].

The overall best scenario is the combination of NBS approaches in both the lower and upper parts of the river. Under the A1B1 scenario, nitrogen export was reduced by 340 tons/year in 2050.

Efforts to reduce phosphorus export are complicated by an expected increase in rainfall intensity in 2050 in relation to the soil characteristics of the Piave river catchment. This increased export is partially mitigated by the A1 and completely reversed in the B1 scenarios. The best scenario is the A1B1 scenario which showed a 148 ton/year reduction in phosphorus export.

To properly estimate the benefits of the proposed NBS, the impacts on other ecosystem services should be considered. Several of these (flood risk reduction, carbon storage and water yield) may be highly significant. However, to explore only those associated with nutrient reduction, we used costs associated with similar NBS projects focused on nutrient reduction to the Venice Lagoon [100]. These projects (

Table 7. The cost of individual NBS projects and the associated annual reduction in N and P in the Veneto region.

Euro	N Tot Ton/Year	P Tot Ton/Year
7,230,000 €	44.9	3.2
4,130,000 €	27.4	1.47
6,600,000 €	43.49	1.4
3,160,000 €	4.9	0.5
825,000 €	1.066	0.45
1,560,000 €	2.14	0.50
2,450,000 €	5.8	0.66
1,650,000 €	4.31	0.66
1,440,000 €	2.56	0.22
1,000,000 €	1.08	0.28
2,066,000 €	18	0.81
1,033,000 €	13.94	1.56
3,430,000 €	5.02	1.08

) had a range of secondary benefits that were not evaluated in the present analysis. The estimated costs for the associated reduction for nitrogen and for phosphorus were calculated using a multiple linear regression analysis (R² = 0.79, p < 0.001).

With an expected intervention cost of EUR 1,128,836 (y-intercept), the cost per ton of N removed by NBS amount was 10.88 EUR/kgN per year. The estimated cost for P was approximately double or 23.83 EUR/kgP per year. These compare well to estimated costs from wastewater treatment upgrades reported by Gratziou and Chrisochoidou [104], who estimated that for nitrogen removal, the cost of a project varies from 4613 EUR/m³ to 488 EUR/m³, with total annual operation cost ranges from 204 EUR/m³ to 17 EUR/m³. The unit cost for P removal in different treatment alternatives ranges from 80 EUR/kgP to 120 EUR/kgP (Bashar et al., 2018). Jabłonska et al. [120] estimated the costs of a hypothetical establishment of wetland buffer zones to reduce the non-point source of N and P to be EUR 9 ± 107 M to remove 11%–82% N and 14%–87% P load from the catchment. This translates into a cost that is saved due to the lack of eutrophication and the problems associated with it [121].

Considering the annual reduction in nitrogen and phosphorus exports achieved using NBS (A1 and B1), compared to scenarios with no NBS (A0, B0), the value of these interventions per year is equal to EUR 2,244,677 for A1 and EUR 2,610,888 for B1, using 2018 values. Given the long-term benefits and multiple risks of climate change, as well as the important secondary benefits, the NBS solutions are well justified.

5. Conclusions

There is growing awareness that NBS can mitigate climate change impacts while securing ESs. The results from the present study, estimating nutrient load and delivery under four different NBS scenarios, suggest that the nutrient load and delivery in the Piave river catchment would benefit from NBS actions, one of which is already underway through rural abandonment. Upper catchment reforestation has increased wooded areas in agricultural lands that are no longer used by local farmers. Reduced water runoff and sediment loss have had and will continue to have positive impacts on sediment, nitrogen and phosphorus loads to the river.

The proposed NBS for the lower catchment (32 km) riverine corridor, based on increased riparian vegetation and erodible areas (A1), will have a more complex activation due to the intense agricultural and urban use that characterizes this part of the catchment. An increase in vegetated areas in the riparian area will reduce the transport of phosphorus-loaded particulates from agricultural soils as a result of erosion. While A1 has a smaller effect on overall nutrient export with respect to B1, this NBS is likely to have more impacts on other major ecosystem services related to biodiversity and recreational value due to the more elevated population density. Decreasing trends in P export by rivers in Europe have mainly resulted from environmental and agricultural policies leading to reduced nutrient inputs to river catchments. On the other hand, nitrate remains a major challenge in many catchments and already compromised receiving waters, such as the Adriatic Sea. The IPCC Climate Change and Land Report [122] emphasizes the need to explore the mitigation potential of restoration actions and the improved management of forests.

It is clear that increased participation in the management and monitoring of our river environments is necessary to improve their status further. Citizen science represented an additional tool to complete the information and for model validation in this study and others. Involving citizen participants directly in monitoring activities can generate a powerful tool to complete the lack of information and improves communities’ influence on management policy in their territory. Understanding the direct and indirect influence of human activities is the first step to evaluate the economic value of nutrients removal. Our estimated cost reflects the gains from investing in NBS and shows how the evaluation of ecosystem services can provide a complete evaluation of where and how watersheds can implement these approaches.