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

: 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.


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.

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

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, Sustainability 2021, 13, 10629 4 of 21 river consortiums and available empirical data. Nutrient flows were divided into sedimentbound (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.

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. 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 −1 y −1 ) 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). 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 3 .

•
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 [ 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 3 -N) and phosphate (PO 4 -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 4 -P concentrations were estimated colourimetrically using inosine enzymatic reactions in seven specific ranges from 0.02 mg L −1 to 1. The complete biophysical table is available in the Supplementary Materials (Table S1).

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:  [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.

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.

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 Sustainability 2021, 13, 10629 7 of 21 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.

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.

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). 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.

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.

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 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-thanusual winter temperatures and changes in the precipitation regimes in much of Europe. 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.   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.   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. 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). 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).

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) was used, which did not differ significantly from measurements made at Nervesa della Battaglia (130 m 3 /s).

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) was used, which did not differ significantly from measurements made at Nervesa della Battaglia (130 m 3 /s). ARPAV samples (Table 3)   Trained citizen scientists collected more than 100 samples along the river in the lower catchment (Table 4), where the number of environmental agency monitoring stations is limited. Nitrate-nitrogen (NO 3 -N) concentrations were high in autumn (1.63 mg/L) and low in summer (0.90 mg/L). Phosphate-phosphorus (PO 4 -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), with surface loads dominating over subsurface loads. T-P 5.96 × 10 2 1.14 × 10 4 -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). 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 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).

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]. 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).

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).
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).

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

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) 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 2 = 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 3 to 488 EUR/m 3 , with total annual operation cost ranges from 204 EUR/m 3 to 17 EUR/m 3 . 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.

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 phosphorusloaded 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.

Data Availability Statement:
The raw data supporting the conclusions of this article will be made available by the authors without undue reservation.