Integrated Valuation of Nature-Based Solutions Using TESSA: Three Floodplain Restoration Studies in the Danube Catchment

Floodplain restoration measures are among the most well-known nature-based solutions for flood risk reduction but practitioners see their limitations in comparison to technical measures when considering both their effectiveness and profitability. The aim of this study is to show the co-benefits (besides flood risk reduction) of floodplain restoration and handle them in terms of monetized ecosystem services (ES). Our work focused on six ES groups for three study areas in the Danube catchment along the Krka, Morava, and Danube rivers. ES mapping through stakeholder engagement is also considered. We applied the methodologies suggested in the Toolkit for Ecosystem Service Site-Based Assessment (TESSA) complemented with alternative methodologies (e.g., questionnaires on social media). Results show annual combined benefits of floodplain restoration in a range from 237,000 USD2019 at Krka to 3.1 million USD2019 at Morava, suggesting the utility of ES assessment. The combination of stakeholder workshops and the TESSA guidelines, as well as the newly developed methods, were all central tools to provide decision-makers with arguments to use nature-based solutions for an integrated and holistic riparian land use management.


Introduction
Floodplains are of great value for humanity and ecology but are worldwide disappearing. In Europe, the World Wide Fund compared the current floodplain areas size with the 19th Century and showed that less than 20% of the former floodplain area remains in the entire Danube River Basin (DRB) [1]. Yet, floodplains cover many functions, such as playing a crucial role in flooding water retention. The International Union for Conservation of Nature (IUCN) defines nature-based solutions (NBS) as "actions to protect, sustainably manage, and restore natural or modified ecosystems, that address societal challenges effectively and adaptively, simultaneously providing human well-being and biodiversity benefits" [2]. In this context, floodplain restoration is an increasingly popular NBS that consists in different soft engineering approaches combining the preservation or improvement of biodiversity and the increase of flood risk resilience. Examples of restoration measures include re-meandering of the main river channel, increasing the floodplain width, and removing artificial dykes or other obstructions to flood flows [3].
Floodplain restoration is often seen as a win-win NBS. For example, in flood risk management, the technical measures used in the last century to protect us against extreme flood events have proved to be not resilient for two reasons. In some cases, the possibility In fact, the fourth step for flood management developing in the direction of integrated approaches is to test innovative flood management strategies in pilot projects [11]. This paper is aimed at filling this research gap, by also testing the same methodology to three pilot areas at the same time.
Our research objective is to develop a systematic method to evaluate ES of restored floodplains based on the well-established TESSA approach. More specifically, we address questions, such as the following: • What are the main benefits of a floodplain restoration scenario? • Does the same method apply for different (but spatially not too far) case studies? • What are alternatives to TESSA's missing methodologies for some ES?
For three cases in the DRB, we tested our approach providing in-depth analysis on relevant ES (carbon storage, greenhouse gases sequestration, flood protection, nutrients retention, cultivated goods, and nature-based recreation) for future riparian land use planning and management. The research provides a new application of TESSA's Version 2.0, especially considering the stakeholder consultation enabling to formulate recommendations for further methodological improvement in times of increasing digitalization and social distancing.

Materials and Methods
To support decision-making processes, we developed a method to assess the potential ES availability of planned but not yet implemented restoration projects. We used the version 2 of TESSA [14] as theoretical background for ES estimation and in some cases for ES valuation. All steps were implemented in a Python code that can be run from QGIS3 [22]. The code and illustrative input data are available on GitHub [23]. For the ES assessment, a consistent quantity of input data was necessary, such as agricultural yields, the population density, or emission factors of greenhouse gases. When lack of data characterized the study areas, we used publicly available data from different institutions and databases (e.g., the Intergovernmental Panel on Climate Change (IPCC) reports [24,25], FAOSTAT [26], EUROSTAT [27,28], and EarthStat [29]). Moreover, the maps resulting from stakeholder workshops were used as input data for the estimation of global climate regulation and cultivated goods ES. In order to evidence potential changes, we applied the methods for both prior restoration state, namely the current scenario (CS) and the restoration scenario (RS).

Study Areas and Restoration Measures
This research is a direct outcome of the Interreg Danube Floodplain project [30] and the methodological approach has been applied to three study sites within the DRB (Figure 1). Begecka Jama (394 hectares) is located at the Serbian Danube, upstream from the City of Novi Sad; the Krka study area (4115 hectares) in Slovenia includes the town of Kostanjevica na Krki, the Krakovski forest, and the Krka river; the Morava study area (17,068 hectares) is located at the border between Czech Republic and Slovakia, and it comprises the area between the Morava and the Thaya rivers, as well as the Morava river itself.
The floodplain restoration measures considered in the RS scenario consist of modifying the water regime to affect the flow conditions and the water supply in the floodplain areas throughout the year. In all study areas, these modifications have the goals of improving flood risk management, protecting the ecosystems, increasing biodiversity, improving habitat quality (terrestrial and aquatic, e.g., spawning areas or better conditions for fish migration), and increasing the diversity and quality of ecosystem services. The floodplain restoration measures considered in the RS scenario consist of modifying the water regime to affect the flow conditions and the water supply in the floodplain areas throughout the year. In all study areas, these modifications have the goals of improving flood risk management, protecting the ecosystems, increasing biodiversity, improving habitat quality (terrestrial and aquatic, e.g., spawning areas or better conditions for fish migration), and increasing the diversity and quality of ecosystem services.

Floodplain Restoration in the Begecka Jama Study Area
Begecka Jama's wetland habitats and hydrological regime have worsened due to forestry activities, pollution from the arable land, flood protection measures, and intensive land use [31]. According to the Danube River Basin Management Plan (DRBMP) update of 2015 [32], this Danube section is a heavily modified water body (HMWB) and its chemical status is "poor", where the standardized scale used within the Water Framework Directive (WFD) of the European Parliament [33] rates as "good" the chemical status of surface waters, whose concentrations of all priority substances do not exceed those permitted by the Environmental Quality Standards Directive 2008/105/EC [34,35]. The area's ecological potential is "moderate" (also derived from the WFD [33], as expressed by biological quality elements, hydromorphological, and physico-chemical parameters [34]). The potential restoration measures analyzed in this study would include a weir reconstruction and diversifying the river morphology and cross profiles (e.g., cleaning and widening the connecting channel between the river and the lake, deepening the existing oxbows and channels, excavating new channels between the deepened oxbows). These measures

Floodplain Restoration in the Begecka Jama Study Area
Begecka Jama's wetland habitats and hydrological regime have worsened due to forestry activities, pollution from the arable land, flood protection measures, and intensive land use [31]. According to the Danube River Basin Management Plan (DRBMP) update of 2015 [32], this Danube section is a heavily modified water body (HMWB) and its chemical status is "poor", where the standardized scale used within the Water Framework Directive (WFD) of the European Parliament [33] rates as "good" the chemical status of surface waters, whose concentrations of all priority substances do not exceed those permitted by the Environmental Quality Standards Directive 2008/105/EC [34,35]. The area's ecological potential is "moderate" (also derived from the WFD [33], as expressed by biological quality elements, hydromorphological, and physico-chemical parameters [34]). The potential restoration measures analyzed in this study would include a weir reconstruction and diversifying the river morphology and cross profiles (e.g., cleaning and widening the connecting channel between the river and the lake, deepening the existing oxbows and channels, excavating new channels between the deepened oxbows). These measures should have the effect of creating new fish spawning areas, improving fish migration conditions and increasing consequently biodiversity, as well as enabling inflow and outflow from the floodplain system. Begecka Jama study area is investigated by the Jaroslav Cerni Water Institute [36].

Floodplain Restoration in the Krka Study Area
The Krka area is protected under the Habitat and Bird Directives, Natura2000, and the International Union for the Conservation of Nature and Natural Resources Red List, also due to the Krakovski forest, which represents the largest lowland floodplain forest in Slovenia [31]. Its overall ecological status was classified as "good" [32]. However, the forest suffers from droughts throughout the year due to the floodplain's disconnection. To solve this problem and improve flood risk management for the Kostanjevica na Krki town, the local authorities [37] hypothesized a combination of measures to increase water conductivity in the river bed, with actions like river bed deepening and the excavation of three corridors to enable floodplain activation [31]. The Krka study area is investigated by the Slovenian Water Agency [37].

Floodplain Restoration in the Morava Study Area
The Morava river faced in the past channel straightening and the interruption of longitudinal continuity by weirs and sills [31]. The river is now an HMWB, with ecological status of class "moderate" and hydro-morphological quality of class "poor" [32]. A restoration scenario was planned, in which certain levees and barriers are removed or adjusted, flood dykes are relocated, river meanders are restored, and oxbows are reconnected to the main Morava channel [38]. The Morava study area is investigated by the Czech Morava River Basin Authority [38] and the Water Research Institute of Slovakia [39].

Stakeholders' Consultation
The stakeholder workshops took place in the study areas between January and February 2019. During the meetings, between 17 (Krka), 37 (Begecka Jama), and 41 (Morava) stakeholders from various interest groups first attended a presentation regarding the state of the floodplains and the potential restoration measures. The stakeholders belonged to different interest groups and went from local to international level. Therefore, it was not possible to differentiate them into groups belonging to the floodplain's upstream and downstream areas. In a second step, the participants were divided into multiple mixed groups and discussed, identified, and marked the locations and intensity of use (on a scale from 1 to 5, where one stands for "Missing to very low" intensity and five stands for "Very high" intensity) of different ecosystems in the current state scenario of the pilot areas. In a third step, the stakeholder groups were asked to repeat the second step, but for a potential floodplain restoration scenario; in this case, the intensity of ecosystem services had to be only recognized as increasing, decreasing, or stable. As a result of the meetings, Table 1 reports the ES that were recognized by the stakeholders in the study areas [40]. A summary of the corresponding potential TESSA methods used for each ES sub-group in this paper can be seen in the last column of Table 1.

Global Climate Regulation
The ES of global climate regulation refers to the storage of carbon and the exchange of carbon dioxide (CO 2 ) and other greenhouse gases (GHGs) between the atmosphere and the plants, the animals and soil within ecosystems. As suggested in TESSA, these were estimated following the Tier 1 methodology of the Intergovernmental Panel on Climate Change (IPPC) reports [24,25] and taking the estimates of rates for CO 2 flux, methane (CH 4 ) flux and nitrous oxide (N 2 O) flux of various habitat types from the above-mentioned sources or from Anderson-Teixeira and DeLucia [42] (Table S1 of the Supplementary Materials).

Carbon Storage
We calculated the total carbon stocks as the sum of the carbon stocks of the aboveground biomass (ABG), the below-ground biomass (BGB), the litter biomass (LB), the dead wood biomass (DWB), and the soil organic carbon (SOC). For Begecka Jama and Krka, we estimated the SOC assuming to have only mineral soil and using IPCC's Tier 1 method. For Morava, due to the presence of organic soil (Luvic Chernozem) as reported in the Food and Agriculture Organization (FAO) Soil Map of the World [43], we extracted the SOC value from a global soil organic carbon map (GSOCmap, v1.5.0) [44]. This makes the process different than TESSA's recommendations.

Greenhouse Gases Flux
We estimated carbon sequestration (i.e., increase of carbon stocks) by assuming that the carbon stock change is taking place only in tree-dominated areas. Differently than what is suggested in TESSA in the Decision tree C6, we calculated the growth of carbon stocks using the growing rates of planted trees, i.e., the Mean Annual Increment (MAI) expressed in annual cubic meters per hectare and taken from the Planted Forests Database (PFDB) [45].
We estimated the losses from the carbon stocks as suggested in TESSA in the Decision tree C7 and in Climate Method M8, i.e., we assumed that carbon losses are derived from disturbances that can come from wood removals, fuelwood collection and charcoal removals, or other disturbances (e.g., illnesses, fires, etc.). The "Forestry Production and Trade" section of the FAOSTAT database [26] provides data on the national level on annual roundwood removals, annual fuelwood removals, and annual charcoal removals in m 3 /year. We then scaled the country-level data from the country values to the study area size. The values used for the reference year 2017 (default year) can be seen in Table S2 of the Supplementary Materials. For Krka and Morava, disturbances (such as illnesses and fires) could not be estimated and were therefore neglected. For Begecka Jama, we used disturbance values from personal communication [46] and assumed that the area in which the disturbance takes place is equal to the whole harvest area.
To estimate the quantity of GHGs emitted/sequestered to/from the atmosphere, we followed the steps suggested in TESSA for CO 2 emissions from drained soils, CH 4 emissions from grazing animals, CH 4 emissions from wetlands (besides for Begecka Jama, for which we had no wetland description), and N 2 O emissions from agriculture. We then calculated their corresponding CO 2 equivalent and overall GHGs flux.

Monetary Value of Carbon Storage and GHGs Flux
We calculated the corresponding monetary value of the stored carbon and of the GHGs flux by multiplying the estimated CO 2 equivalents times the values of the CO 2 emissions taxation systems documented in the report of the World Bank [47]. The Slovenian Carbon tax rounded up to the nearest integer is 19 USD 2020 per metric tons of carbon dioxide equivalent (tCO2e) [47], as well as the European Union (EU) Emissions Trading System (ETS), for year 2020 [48]. In the previous years, the EU ETS values were 16 USD 2020 per tCO2e in 2018 and 25 USD 2020 per tCO2e in 2019 [48]. Since the overarching framework of the international carbon market remains unclear and decisions for future prices in the EU are postponed to 2021 [47], we used the values from 2018 and 2019 to estimate error calculations of the values of stored carbon and GHGs flux services.

Water-Related Services: Flood Protection
The reduction of flood risk was originally one of the main goals of the floodplain restoration NBS and of the Danube Floodplain project [30]. We therefore put particular attention in the estimation of this ES. First, we collected three water depth maps for each study area, resulting from the 2D hydrodynamic modeling of the three return period (RP) groups of high probability (RP 2 to 5 years), medium probability (RP 10 to 20 years), and low probability (RP 100 years). The modeling was conducted under the framework of the Danube Floodplain project [30] for both CS and RS scenarios. Secondly, we applied the Joint Research Centre (JRC) damage functions [41,49] to all scenarios and land uses to estimate the flood-caused damage in our study areas. The damage curves can be observed in the Supplementary Materials. Finally, we applied the trapezoidal method for flood risk (expected annual damage, EAD) estimation [50], as shown in the following function: where n = 3 is the number of return periods, T is the return period in years (shown in detail for each study area in the Supplementary Materials, together with their corresponding lower and upper uncertainty boundaries), and D is the corresponding damage.

Water Related Services: Nutrients Retention
Although some steps overlap with the guidelines, the estimation of the nutrients retention by the floodplains did not follow TESSA because we did not have access to measured data of water quality upstream and downstream of the studied floodplain areas. Instead, we analyzed the data from the DanubeGIS [51] of total nitrogen (TN) measurements at the Danube and its tributaries and combined them with our knowledge on the presence of active floodplains in the DRB [30]. We analyzed comparable measurements To understand the TN retention of the whole floodplain, i.e., to scale the value in mg N/l/ha to a total value of retained kg TN, we needed the volume of water filtered by the floodplain per year. Therefore, we took the floodplains' activated volume that we simulated for extreme flood events (HQ2 to HQ5, HQ10 to HQ20, and HQ100) and calculated the expected annual retention volume (EARV) with the trapezoid method, as shown in the following function: where = 3 is the number of return periods, is the return period in years, and is the corresponding retention volume. The specific values of (together with their corresponding lower and upper uncertainty boundaries) and for each study area can be found in the Supplementary Materials. The estimation is valid under the assumption that the volume that is additionally retained by the restored floodplain in comparison to the CS scenario is also the volume that is additionally filtered by the floodplain.
To attribute a monetary value to the TN retention of the floodplain, we applied the benefit transfer (BT) method by using the database of floodplains' ES values in the DRB and its intersecting countries [9]. We used only the values expressed in USD 2018 /kg N and applied their average of 7.27 USD 2018 /kg N ( Figure 2) to the estimated annual quantity of retained TN for each study area. To estimate the corresponding errors, we applied the values 3.75 USD 2018 /kg N and 10.79 USD 2018 /kg N, being these the first and third quartiles of the benefit transfer values, respectively.

Cultivated Goods
The cultivated goods provisioning ES includes three types of goods: crops, livestock, and aquaculture. These ES were calculated using a combination of interviews to local authorities and publicly available data, such as FAOSTAT tables for agriculture and livestock [26] and EUROSTAT tables for aquaculture [28]. The crops goods estimation additionally required the EarthStat maps of harvested areas and yield for each crop type [29]. With this information, it was possible to calculate the annual total yield of each listed crop type for the selected areas. The knowledge on the crops, livestock, and fish species present in the study areas was provided by the interviews [46,52,53] and is shown in Table S7 of the Supplementary Materials. The spatial extension of the production areas was given instead by the stakeholder workshops ES maps. The national data were scaled according to the size of the area. The cultivated goods ES values were estimated with market prices. The data used were statistics from 2017 of the "Trade-Crops and livestock products" section of the FAOSTAT database [26], which provides the producer prices per unit, or the To understand the TN retention of the whole floodplain, i.e., to scale the value in mg N/l/ha to a total value of retained kg TN, we needed the volume of water filtered by the floodplain per year. Therefore, we took the floodplains' activated volume that we simulated for extreme flood events (HQ2 to HQ5, HQ10 to HQ20, and HQ100) and calculated the expected annual retention volume (EARV) with the trapezoid method, as shown in the following function: where n = 3 is the number of return periods, T is the return period in years, and RV is the corresponding retention volume. The specific values of T (together with their corresponding lower and upper uncertainty boundaries) and RV for each study area can be found in the Supplementary Materials. The estimation is valid under the assumption that the volume that is additionally retained by the restored floodplain in comparison to the CS scenario is also the volume that is additionally filtered by the floodplain.
To attribute a monetary value to the TN retention of the floodplain, we applied the benefit transfer (BT) method by using the database of floodplains' ES values in the DRB and its intersecting countries [9]. We used only the values expressed in USD 2018 /kg N and applied their average of 7.27 USD 2018 /kg N ( Figure 2) to the estimated annual quantity of retained TN for each study area. To estimate the corresponding errors, we applied the values 3.75 USD 2018 /kg N and 10.79 USD 2018 /kg N, being these the first and third quartiles of the benefit transfer values, respectively.

Cultivated Goods
The cultivated goods provisioning ES includes three types of goods: crops, livestock, and aquaculture. These ES were calculated using a combination of interviews to local authorities and publicly available data, such as FAOSTAT tables for agriculture and livestock [26] and EUROSTAT tables for aquaculture [28]. The crops goods estimation additionally required the EarthStat maps of harvested areas and yield for each crop type [29]. With this information, it was possible to calculate the annual total yield of each listed crop type for the selected areas. The knowledge on the crops, livestock, and fish species present in the study areas was provided by the interviews [46,52,53] and is shown in Table S7 of the Supplementary Materials. The spatial extension of the production areas was given instead by the stakeholder workshops ES maps. The national data were scaled according to the size of the area. The cultivated goods ES values were estimated with Sustainability 2021, 13, 1482 9 of 22 market prices. The data used were statistics from 2017 of the "Trade-Crops and livestock products" section of the FAOSTAT database [26], which provides the producer prices per unit, or the EUROSTAT, which provides the prices at first sale for human consumption in EUR/year [28]. To estimate the results' uncertainty boundaries, we used the minimum and maximum national statistics values of primary production (for livestock goods), producer prices (for agricultural goods) or both (for aquaculture goods) in the periods 2014 to 2018 (for agricultural and livestock goods) or 2008 to 2017 (for aquaculture goods).

Nature-Based Recreation
To assess the nature-based recreation (e.g., exercising, experiencing nature, etc.) provided by the floodplain areas and their restoration, we applied the individual travel cost method (ITCM). As a response to the COVID-19 pandemic and its consequent travel restrictions [54][55][56], this method was based on interviews that were conducted online through LimeSurvey [57] from 7 August 2020 to 1 September 2020. We used Facebook events [58] and Instagram [59] posts (with hashtags related to the study areas) to advertise the survey (in locations with a radius of 20 km around Begecka Jama, 20 km around Kostanjevica na Krki, and 40 km around Lanzhot for Morava). To retrieve data on the restoration scenario, the interviews included a section in which the respondents described their potential reaction to a hypothetical floodplain restoration. A template of the interviews can be found in the Supplementary Materials. The ITCM requires as input data the count of the visits of an individual to a site in a year, the corresponding travel cost (TC) to the site (sum of cost to get to the site with fuel prices for each country from the European Commission [60] and additional expenses), and can include other characteristics (e.g., age, education level, etc.).
where α is the intercept, and β and γ are the coefficients estimates. Based on the fitted Poisson model, the consumer surplus per visit was calculated as the negative inverse of the constant (−1/β) of the TC variable. Herein, the consumer surplus is defined as the difference between the total travel costs sustained by a visitor to a recreational site and the maximum amount the visitor is willing to pay to make a visit to that site. Multiplying the consumer surplus by the total number of visits gave a total consumer surplus for the site [63]. To estimate the results' uncertainty boundaries, we propagated the lower and upper boundaries derived from the standard error of the β coefficients. The total number of visits was retrieved from additional e-mail conducted interviews [64][65][66][67] and personal communication with local authorities [53].

Stakeholders' Consultation
The combination of the stakeholder workshops and of an expert-based analysis of land uses [68,69] resulted in six digital maps, one for the current state scenario and one for the restoration scenario in three study areas. These are in shapefile format and include in their attribute tables the corresponding ES types and intensities. The maps are shown in the Supplementary Materials.

Global Climate Regulation
The results of carbon storage are presented in the first part of Table 2. For all three study sites and for both scenarios, the largest carbon stock is represented by soil organic matter and the smallest by either BGB or litter and dead wood carbon. As a result of forest-or grass-dominated area conversion as alternative to crop-dominated areas, carbon storage would increase by 12.5%, 1.8%, and 2.9% for the Begecka Jama, Krka, and Morava areas, respectively. This would lead to a win of stock monetary value in the area of 0.4 (for Begecka Jama) to 4.0 (for Morava) million USD 2020 as a total static value (not annual). The second part of Table 2 shows the net GHGs fluxes in the current and restoration scenarios, where the negative values represent equivalent CO 2 emissions and positive values represent equivalent CO 2 sequestration. The results are dominated by CO 2 emissions due to the carbon stock losses, if the tree-dominated areas would be considered to be harvested or affected by other disturbances. Emissions of N 2 O and CH 4 also show substantial effects on the overall GHGs balance in both states. In the restoration compared to the current scenarios of Begecka Jama, the corresponding increase of sequestration of equivalent CO 2 is 14.9%. On the contrary, the additional presence of wetlands in Krka causes an increase in CO 2 emissions and decreases the presence of forests, which consequently decreases the carbon stock increment in RS, resulting in the end with a negative effect of the RS on equivalent CO 2 emissions with an increase by 2.1%. In the Morava area, the sequestration of equivalent CO 2 increases by 8.0% for the restoration scenario due to an increase of carbon stock and a decrease of methane emissions.

Water Related Services: Flood Protection
In terms of flood protection, the NBS tested in the three study areas show inhomogeneous results. As Table 3 shows, at Begecka Jama the annual value of flood protection provided by the NBS within the 2D model zone was estimated to decrease by 0.07%, a value that could have also been affected by the accuracy of the estimation itself. The modelled flood risk benefits of flood storage resulting from floodplain restorations in Krka and Morava are instead positive (respectively, by 2.16% and 43.67%), though merely affected by the small difference in water level. The floodplain's ES flood protection is shown in terms of the expected annual flood-caused damage, i.e., not in terms of ES benefits.

Water Related Services: Nutrients Retention
The results on ES value of total nitrogen (TN) retention (Table 4) depend on the amount of filtered volume and on the flooded area size. For Begecka Jama and Krka, where either the filtered volume or the flooded area of tree-, grass-, and wetland-dominated areas is increasing after implementing the potential NBS, the ES value is also increasing (by 6.2% and 9.7%, respectively), while the flooded area is decreasing in Morava with the NBS, and its corresponding ES value decreases by 8.0%.

Cultivated Goods
The results of cultivated goods' provisioning ES value for crops, livestock, and aquaculture are reported in Table 5. Both scenarios of Begecka Jama include no cropland nor aquaculture, while the livestock products provisioning value would decrease by 2668 USD 2017 /yr, in case of the floodplain restoration measures. In Krka, we estimated that the mean annual net benefit from this type of provisioning ES would decrease by 5384 USD 2017 /yr, due to the substitution of a part of the area's cropland and grasslands with wetlands. In Morava, the tree-dominated areas would increase and would cause a lower presence of livestock with its consequent lower revenue (by 458,309 USD 2017 /yr) for the study area's cultivated goods ES value.

Nature-Based Recreation
The precision level was calculated as suggested by TESSA [14], by considering the results of the first ten surveys for each study area. The resulting precision level is approximately 80% for Begecka Jama, 65% for Krka, and 55% for Morava. This corresponds to sample sizes of 134 for Begecka Jama, 47 for Krka, and 50 for Morava. Based on the interviews for the census of visitors to the area, we used as number of visitors per year to each study area of 10,000 for Begecka Jama (only given value), 25,000 for Krka (average of two values), and 100,000 for Morava (median of three values). Table 6 shows the corresponding results. All three study areas show an increase of the recreational ES value in case the NBS would take place because the visit numbers would increase, since we kept the consumer surplus the same for CS and RS. The benefits from recreation would increase by 102%, 5%, and 20% for Begecka Jama, Krka, and Morava, respectively. A representation of the respondents' age is shown in Figure 3, where we can see all generations represented. values), and 100,000 for Morava (median of three values). Table  ing results. All three study areas show an increase of the recreat NBS would take place because the visit numbers would increas sumer surplus the same for CS and RS. The benefits from recr 102%, 5%, and 20% for Begecka Jama, Krka, and Morava, respe of the respondents' age is shown in Figure 3, where we can see all

Sum of All ES
The total absolute value of the ES benefits (i.e., the sum of R of the three floodplain restoration measures (without conside stocks) was estimated at approximately 1.5 million USD 2019 /yr USD 2019 /yr in Krka, and 3.1 million USD 2019 /yr in Morava. This carbon storage values because these are not expressed in mon only in static monetary value.
However, the total added value of the NBS in terms of ES is

Sum of All ES
The total absolute value of the ES benefits (i.e., the sum of RS-CS in the above Tables) of the three floodplain restoration measures (without considering the gains in carbon stocks) was estimated at approximately 1. the most affecting ES types (flood protection, GHGs flux regulation, cultivated goods, and nature-based recreation) are not constant among study area, nutrients retention is constantly the least contributing ES to the sum of benefits (i.e., the sum of RS-CS in the above Tables).    The bar plot in Figure 5d also displays the NBS added value per unit area. Krka and Morava show comparable trends, in which the added ES value are of the same order of magnitude for flood protection (1) and nature-based recreation (2). The latter has instead the highest value per area unit for Begecka Jama, which is also mainly profiting from the RS in terms of GHGs sequestration and nutrients retention. The bar plot in Figure 5d also displays the NBS added value per unit area. Krka and Morava show comparable trends, in which the added ES value are of the same order of magnitude for flood protection (1) and nature-based recreation (2). The latter has instead the highest value per area unit for Begecka Jama, which is also mainly profiting from the RS in terms of GHGs sequestration and nutrients retention.

Analysis of Singular ES Groups
In terms of carbon storage, all three study areas benefit from the restoration measures, mainly due to increase of carbon biomass above ground and from organic soil. We observed that the larger the study area's size, the higher the absolute value of this ES. However, we consider carbon storage ES value as an expression of the natural capital, and not as an ES, from which beneficiaries can get annual revenue from the floodplain. Therefore, we did not include this value in the total added ES value of the NBS.

Analysis of Singular ES Groups
In terms of carbon storage, all three study areas benefit from the restoration measures, mainly due to increase of carbon biomass above ground and from organic soil. We observed that the larger the study area's size, the higher the absolute value of this ES. However, we consider carbon storage ES value as an expression of the natural capital, and not as an ES, from which beneficiaries can get annual revenue from the floodplain. Therefore, we did not include this value in the total added ES value of the NBS.
In terms of GHGs fluxes, results are heterogeneous among the study areas. The estimation of the disturbances in tree-dominated areas in Krka has a strong negative effect on the total. In general, the assumptions made for carbon stock losses have a great influence on the total GHGs balance. The uncertainty in tree disturbances are added to the ones of the IPCC Tier 1 method, which is a standard and very useful tool but which is also highly sensitive to the input data chosen. Therefore, in case time and resources would be available, we would have chosen a more accurate method (e.g., including field measurements), as suggested by TESSA.
The value of the NBS in terms of flood risk protection is noticeable for all three pilot areas mainly in the absolute value. The results are comparable to the results of Peh et al. (2014) [19], who found the net flood protection of a wetland restoration of 23,075 USD/yr. However, when looking at the avoided flood risk maps, these differences are low (if not even negligible) because the difference in water level (used as input data) is smaller than one meter. In addition, the flood risk estimation is highly affected by the stepwise damage function chosen for the estimation and by which land uses are recognized as damage-prone by the damage function type itself (e.g., the Begecka Jama study area did not include any residential land uses). As a result, the higher value for Morava is shown to be a factor of the size of the study area, more than the effectiveness of the measure itself.
The output of cultivated goods added ES value are always negative, as expected since in all three study areas the area available for agriculture or livestock would decrease in in case of a floodplain restoration. In comparison to other studies [18,19,21,[70][71][72][73], the losses of cultivated goods values due to land use change are much smaller because the decrease in agriculture or livestock suitable area is also quite low. This is an encouraging output of the presented NBS because stakeholders, especially in European countries, often associate floodplain with high losses in terms of agricultural production.
The accuracy on the estimations of cultivated goods depends on the quality of the data on crop land area and yields provided by the EarthStat dataset [29]. This has a cell grid resolution of 0.0833 degrees, corresponding to a minimum two and maximum nine cell grids. The availability of local data, e.g., through interviews, could probably increase the accuracy of this estimation. However, for a rapid assessment for decision-makers, we consider sufficient the mixed approach used in this study. In fact, we used input from the stakeholders to derive the most important crop types and input from publicly available sources, to get the market prices, which are country-specific.
Regarding nature-based recreation services, we found the results on the consumer surplus to be reliable and comparable to previous studies of nature-based recreation [18,19,21,63,70,[72][73][74], where the net value of restoration was found to be between around 60,000 USD/yr [73] and more than 4 million USD/yr. In our case, the highest value was found for the Krka area, which could be explained by the high range of nature-based activities offered in the area and including the city Kostanjevica na Krki, a main Slovenian tourist attraction close to the capital city Ljubljana. Similarly, Begecka Jama is located in a strategic position, reachable in less than two hours' drive from two big cities (Novi Sad and Belgrade), and its surplus value is almost as high as Krka's. The lowest value of the consumer surplus in Morava might have been affected by the presence of so many other touristic attractions in proximity that could compete with each other. However, these values were calculated without considering the potential increase of mosquitos' presence, which could reduce the recreational amusement in the summer months. Additionally, our study assumed that the areas would remain easily accessible after the restoration measure would be implemented.
The final ES value of nature-based recreation is strongly influenced by the number of visits that take place annually, a value which is difficult to estimate, since areas in proximity are also natural areas. Nevertheless, the interviewed people showed a higher visitation rate, in case the NBS would be implemented. In addition, to increase the precision level, the number of interviews might be increased, although, according to some methods (e.g., published tables in Israel [75]), even for a population size higher than 100,000, the sample size for 10% precision would be still 100. We consider the online interviews valuable for our purposes. In fact, the age range of the responses shown in the histogram of Figure 3 proves that the use of online surveys did not necessarily exclude older generations because Facebook and Instagram are not anymore only an instrument of younger generations. The Pew Research Center (2019) [76] showed that U.S. Boomers (born 1946-64) and Silents (born 1945 and earlier) have both increased their Facebook use between 2015 and 2019 from 43% to 60% and from 22% to 37%, respectively. However, we are aware that real-person interviews might have given different results.

Main Outputs
Taking into account six ES, floodplain restoration is bringing added monetary value for all three study areas that we considered. A surprising result is that the total increase of ES of the Krka study area is much smaller in comparison to Begecka Jama and Morava. This can be explained by the fact that the used study area had already a much higher ES level in the current situation.
It is also interesting to note that, as expected, the restoration projects have a different impact to different types of ecosystem services. The provisioning ES (here represented by the cultivated goods) are decreasing in all study areas, while the regulating and cultural services are increasing in a much more complex spectrum of services. These results are in line with previous results from floodplain restoration analyses in Nepal by Merriman et al. [18] and in the U.K. by Peh et al. [19]. As an initial approach for ES-based evaluation of NBS at the Danube and its tributaries, the results can be the basis for a further analysis of the interaction among ES, such as the nexus analysis approach suggested by Fürst et al. [77] and Babí Almenar et al. [78]. This could help us better understand the cause-effect relationship of benefitting from one ES group (e.g., provisioning) to the availability of other ES groups (e.g., regulating or cultural).
The information acquired during the stakeholder workshops was of great help to collect the necessary information about the areas and to map the ES. The results of these workshops were used as input data for the ES assessment with TESSA's methodologies. Although we did not specifically differentiate between the floodplain area's upstream or downstream stakeholders (who would, e.g., benefit from the flow of water services), we were able to bring together different views because the stakeholders' representation was a mixture of local and national public authorities, sectoral agencies, NGOs, and international organizations, and in the case of Krka even the general public. Nevertheless, a broader consultation may have described and judged the ES differently [18].
Looking at the effect of the singular ES types into the sum of benefits, the most affecting ES type varies according to the study area. In Begecka Jama and in Krka, the most affecting ES is recreation, namely the nature-based recreational activities, while in Morava, GHGs flux sequestration has the highest weight, followed by nature-based recreation benefits. This behavior is probably due to the fact that, on the one hand, GHGs sequestration was estimated based on the surface area, which is the highest for Morava. Nature-based recreation ES values are site-specific and depend on social behaviors rather than the size of the study area and are therefore not necessarily proportional. On the other hand, the nutrients retention services showed the lowest effect of this NBS's benefits in all study areas. This could be affected by the ES value transfer or by the approximate method used to estimate the amount of retained nitrogen. To our knowledge, no other publications reported the results on nutrients retention with TESSA. In our case, the floodplain restoration NBS cannot be justified for merely flood risk purposes, especially in Begecka Jama, where the flood protection service is slightly decreasing due to the restoration measure. However, flood protection still brings Krka's second and Morava's third highest benefit contribution to the total services' improvement.

Quality of the Methodological Approach
Geographic information capacity plays a significant role in understanding ES processes [79] and in finding the potential ES hotspots and lowspots of restoration projects. Therefore, general actions to improve the ES assessment at the local level might involve creating a standardized GIS version of the TESSA models, to represent its results spatially. These could be refined for specific regions, e.g., by using local community knowledge. We considered the implementation of the TESSA methodology on a python script written for QGIS as crucial. Once the script was finally written, this choice allowed including input data from freely available sources, but it also decreased the execution time of TESSA tasks. Our technique shows a clear advantage both over the mere mapping of the floodplain's services and over other more time-demanding (e.g., InVEST) software.
TESSA is a helpful tool. The guidelines gave a clear overview of the necessary steps to follow for a quick ES estimation in the study areas. Although the steps are clear and easily implementable, the collection of the big amount of input data is highly time intensive and requires many resources and contacts to local authorities. On the one hand, we encourage the utilization of TESSA for the further evaluation of other kinds of NBS. On the other hand, we invite TESSA's developers to complement methodologies of ES assessment, e.g., by adding guidelines for online interviews or by adding the possibility of using social media, not only for data collection, but also for the design of NBS.
Although we consider the results of TESSA's application useful for a preliminary evaluation of NBS, we found some points of potential improvement. Firstly, in our results, we show the need for TESSA to add more ES within the tool, specifically concerning habitat services.
Secondly, as also recognized by Merriman et al. (2018) [18], the nature of the tool makes TESSA prone to represent mainly those ES that are the easiest to monetarize. Without straightforward methods, the other ES are therefore in danger of being overlooked and under-represented. Accordingly, our results include all ES for which we found readily available methods to estimate their monetary values because we wanted to use a common unit of measure for comparing the scenario and the study areas. This means that we neglected other ES, for which no available methods or data existed, such as noise regulation or local climate regulation. We also encountered difficulties in estimating harvested wild goods, due to a high demand of data. We decided to exclude these ES from our estimation, also reinforced by the assumption that the floodplain restoration would have a low impact on the mentioned ES. In fact, noise regulation and local climate regulation are two ES, which would most likely be affected in urban areas but not by much in our rural study sites. However, stakeholders recognized noise regulation as a floodplain service during the workshop in the Morava study area. In addition, the change in the amount of harvestable goods is very unpredictable for the relatively small changes in our study areas, proven by the fact that stakeholders had heterogeneous and weak opinions on the consequences of the RS with regards to harvested wild goods.
Moreover, we agree with Merriman et al. (2018) [18], who judged the methods suggested by TESSA to assess water quality as too coarse or too time consuming. Nevertheless, we do not want to undermine the importance of investing time and resources in a proper estimation of ES, and we recognize that sometimes an easy and quick solution is not possible to understand complex phenomena.
ES values corresponding to nutrient retention have the lowest effect on the total ES valuation. The methodology used is a new suggestion for the TESSA toolkit, in case no available measurement data and no modelling resources would be available. Data from other studies could also be used as source of information. For example, Doll et al. (2020) [80] found out that for their urban stream restoration project, on average 9% to 15% of the total annual streamflow volume accessed the floodplain, but the percentage of annual streamflow volume that was potentially treated ranged from 1.0% to 5.1%.
As for other studies conducted with TESSA [19], the missing ES quantitative estimations lead to a more conservative result. On one hand, the inclusion of stakeholders in the estimation process allows to include the qualitatively indicated ES in the decision-making process. On the other hand, a bigger picture including all non-monetized ES would be preferable. As an example, the added value of cultural ES and of pollination ES were not included in the estimation, due to difficulties in monetarization for the former, and challenges in knowing about pollinators in the areas for the latter. These factors could potentially be included by a higher engagement of stakeholders, as done by Pugliese et al. (2020) [6].
We also want to underline that the scale of the estimation is highly affecting the accuracy of the results. In contrast to other river-related disciplines (e.g., hydrological modeling, hydrodynamic modeling), the estimation of ES at the local scale is made more difficult the smaller the study area gets and remains a task with high complexity. The biggest difficulties were encountered by the data collection. For example, the application of the FAOSTAT data at the national level would be more appropriate for catchment scales, other than for floodplain scales. In addition, in another example, i.e., flood-caused damage estimation, national-level data were used in form of the flood-damage functions, which could have been more accurate, if local damage or exposure data would have been available.
In this respect, it should be considered that our findings are based on a limited amount of local-specific data.
Besides the above-mentioned phenomena, several other steps of our work are affected by uncertainty, such as the application of the benefit-transfer function, the fit of the Poisson distribution to estimate the visitation rates as a function of the travel costs, and the timeframe used for the monetary values. Herein, we presented a first attempt of error estimation. However, dealing with many variables for different ES, there are even more input variables that should be considered to provide a meaningful error estimation. To fulfill this task, we should to put up a new system to consider all possible sources of uncertainty of the ES estimations, e.g., by using a Monte Carlo simulation. Moreover, most of this uncertainty does not affect the overall results which present the percentage change for each ecosystem service between the two states. For each metric, the error should be similar for both the current state (CS) and restoration state (RS) [21].
Moreover, to describe the uncertainty associated with our ES estimates, we followed TESSA's recommendations and categorized the confidence of the results, choosing among "high", "medium" or "low". Based on these standards given by TESSA, we rated the estimation confidence level of flood protection ES and nature-based recreation ES as "high", that of carbon storage and the GHGs flux services as "moderate", and the confidence level of cultivated goods and nutrients retention services as "low". Therefore, decision-makers should remember the presence of uncertainty, when using these results for decision-making.

Conclusions
We estimated the benefits of floodplain restoration in terms of monetized ecosystem services (ES) in three study areas of the Danube catchment, by also including stakeholder engagement. The conclusions of the research are threefold. We estimated the added value of the benefits of river and floodplain restoration to test the quality and effectiveness of the NBS. We showed that the planning of NBS for flood risk management should not only use standard methods (e.g., hydrodynamic modeling) to support decision-making, but also assess ES for a more holistic picture of the potential consequences of the potential NBS. We provided an example with a mixed application of TESSA and alternative methods, which could be considered by TESSA developers to be included in a new version of the toolkit.
We estimated a total gain of ES in the three pilot areas of approximately 1.5 million USD 2019 /yr in Begecka Jama, 237,000 USD 2019 /yr in Krka, and 3.1 million USD 2019 /yr in Morava. The results are mainly affected by GHG fluxes and changes in nature-based recreation. However, we have shown that the tested NBS will only weakly affect the retention of nutrients. We also observed a diverse effect of the NBS on flood protection among the case studies. Although we did not find that floodplain restoration NBS can be justified for flood protection only, we recall that this output is only valid for these specific study areas and that NBS remain a flexible and resilient way to address natural hazards [81,82].
In this way, we brought further evidence in favor of floodplain restoration measures to be implemented for a general benefit of the communities. In fact, without considering the benefits of NBS, floodplain restoration measures would have much lower chances of being accepted by decision-makers and stakeholders. ES assessment can be useful for decision-makers to locate where to build or restore ecosystems [79]. Policymakers and researchers should give stakeholders a greater role in the design of floodplain restoration measures and in their evaluation, including ES assessment and monetarization. At the same time, researchers should develop new methodologies to rapidly evaluate the missing ES types, which are not included in commonly used ES assessment guidelines (TESSA) or software (InVEST, ARIES, etc.). Moreover, scientists should study the effects of upscaling local-scale methods to the national (or river basin) extent, especially in case more floodplain restoration measures would be implemented at the same time (e.g., nature-based recreation). We suggest that monitoring is done, in case restoration measures would be implemented, to confirm or discard the ES assessment's results. In fact, a high number of factors influence the floodplain ecosystem and the phenomena taking place in it. Therefore, we should make sure that the assumptions used do not invalidate the ES assessments.
Although some progress has been made using our methodology, this approach, being based on a toolkit for rapid execution, assesses only a part of all ES potentially provided by floodplains. In addition, an improvement of the interpretation of the results might be given by analyzing the results' uncertainties. Further, more modeling could be implemented to get a more detailed estimation of some ES (e.g., of water quality).
We finally call for a better inclusion of ES assessment in the Danube River Basin Management Plans, for not only improving ES themselves but also because ES improvement intersects with the achievement and monitoring of the Sustainable Development Goals. ES assessment would act for different purposes, such as to encourage a sustained, inclusive, and sustainable economic growth (Goal 8) and to facilitate a sustainable management of water (Goal 6) and terrestrial ecosystems (Goal 15) [83].
Supplementary Materials: The following are available online at https://www.mdpi.com/2071-105 0/13/3/1482/s1, Section S1: Ecosystem services mapping, Section S2: Additions to the methodology. Data Availability Statement: Publicly available datasets were analyzed, as cited throughout the study. The data resulting from personal communication and questionnaires (see acknowledgements), and the input data used for flood risk estimation (i.e., two-dimensional hydrodynamic modeling results) are not publicly available. Other data presented in this study are available upon request by contacting the corresponding author.