On the Management of Nature-Based Solutions in Open-Air Laboratories: New Insights and Future Perspectives

Abstract

The adoption of Nature-Based Solutions (NBSs) represents a novel means to mitigate natural hazards. In the framework of the OPERANDUM project, this study introduces a methodology to assess the efficiency of the NBSs and a series of Open-Air Laboratories (OALs) regarded as a proof-of-concept for the wider uptake of NBSs. The OALs are located in Finland, Greece, UK, Italy, and Ireland. The methodology is based on a wide modeling activity, incorporated in the context of future climate scenarios. Herein, we present a series of models’ chains able to estimate the efficiency of the NBSs. While the presented models are mainly well-established, their coupling represents a first fundamental step in the study of the long-term efficacy and impact of the NBSs. In the selected sites, NBSs are utilized to cope with distinct natural hazards: floods, droughts, landslides, salt intrusion, and nutrient and sediment loading. The study of the efficacy of NBSs to mitigate these hazards belongs to a series of works devoted to the implementation of NBSs for environmental purposes. Our findings prove that land management plays a crucial role in the process. Specifically, the selected NBSs include intensive forestry; the conversion of urban areas to grassland; dunes; marine seagrass; water retention ponds; live cribwalls; and high-density plantations of woody vegetation and deep-rooted herbaceous vegetation. The management of natural resources should eventually consider the effect of NBSs on urban and rural areas, as their employment is becoming widespread.

1. Introduction

Nature-Based Solutions (NBSs) to societal challenges consist of solutions deeply inspired and supported by nature [1,2]. This concept has been exploited in climate mitigation and adaptation challenges, but it has been recently extended to other fields such as cities [3,4,5,6], land usage [7], and soft science [8]. In the framework of natural hazards mitigation policies, NBSs are receiving increasing attention [9,10,11,12,13,14]. In this perspective, the latest research works are also focusing on specific NBSs effects (e.g., [15,16,17,18]) with the aim of delivering general guidelines for their implementation and management. Eventually, other recent works are focusing both on the institutionalization [19] and the implementation [20] of NBSs to identify their formal description in different frameworks.

The strong interdisciplinarity surrounding the science of NBSs [21] reflects the complexity of their purpose. This is particularly true in the climate change perspective, where the research for solutions must collide with the social, economic, and political aspects simultaneously (e.g., [22,23]). Several studies have been conducted to investigate the effectiveness of the NBSs accounting for the co-benefits for the stakeholders and the resulting trade-off reduction [24,25]. Moreover, NBSs are increasingly seen as a solution for the achievement of several sustainable development goals by bringing a set of societal benefits [26]. These aspects make the NBSs an innovative and challenging set of tools.

In this context, the Horizon 2020 OPERANDUM project aims to reduce hydro-meteorological risks in European territories through green, blue, and hybrid types of NBSs. It aims at the provision of science-evidence for the usability of NBSs; best practices for their design based on participatory processes; it foresees multiple levels of stakeholders’ engagement acceptance to promote its diffusion as a good practice; and it establishes the framework for the strengthening of NBS-based policies. To these aims, the OPERANDUM project delivers methods for the validations of NBSs to mitigate hydro-meteorological hazards across European rural and natural territories, employing several innovative study areas (Open-Air Laboratories—OALs). The concept of OALs has been first introduced in a specific project devoted to the delivery of science products to the community to investigate environmental issues [27,28]. Eventually, a hydrological OAL with similar OPERANDUM-OALs features has been established in Austria in the past years [29,30]. Herein, OALs allow a systematic analysis of the most appropriate mitigation policy, also in terms of cost–benefit analysis. They must not be intended as merely laboratory sites, but locations where real natural hazards are coped with efficient tools. This is possible through the inclusion in the designing and managing stages of a multitude of stakeholders. OALs represent the bond between the decision-makers and the scientific community.

OALs are properly selected around Europe where the green and blue and hybrid types of NBS are installed or are likely to be installed. A series of experiments and measurement campaigns in these locations, different numerical techniques, and interactions with local socio-economic aspects will allow the development, management, and monitoring of selected NBSs depending on the predefined hydro-meteorological hazards (floods, droughts, coastal erosion, nutrient loading, landslides, etc.).

Thus, NBSs belong to the class of natural hazard coping strategies. Several approaches have been adopted to face this issue, including the implementation of stationary processes [31]; remote sensing tools [32,33]; social vulnerability mapping [34,35]; and cost analysis [36]. NBSs represent an innovative tool that will provide additional information to the latter strategies. Indeed, the advance here is represented by the broadening of the local space-scale of the NBSs to the general natural hazard mitigation policies. This is crucial since the policymakers play a fundamental role in the actuation of mitigation policies (see the case of the Rosia Montana mine, reference [37,38]). In this view, NBSs can be seen as helpful tools to policymakers for the understanding of the best actions to undertake to face a specific natural hazard.

Within this frame of reference, a leading role is played by the study of the long-term efficiency of NBSs. It generally depends on a multitude of aspects at different time and space scales that must be simultaneously considered. The lack of information and data at different spatial and temporal scales undermines the assessment of their multidimensional impacts and not least in the future climate condition. Although the NBSs are small compared to the areas of interest, the assessment of their efficiency through numerical models requires a large series of data, ranging from local to regional scale at different spatial and temporal resolutions. Hence, different resolution models coupling and downscaling techniques are needed [14]. This task can be achieved with sequences and/or schemes of numerical models, following a procedure similarly applied in the context of climate modeling (e.g., [39,40]). The proper coupling and merging of the numerical models represent a crucial role because of their multidisciplinarity, applicability to different domains and resolutions, and interoperability. Studies on these aspects are lacking, probably due to the recent development of the NBS concept and application.

In this work, we will present several models’ chains for the estimation of specific NBSs’ efficiency in distinct OALs. The selected OALs are located in five different European countries (Greece, Ireland, Finland, UK, and Italy), which are mainly exposed to floods, droughts, eutrophication, shallow landslides, and erosion hazards. In UK and Finland OALs, the NBSs have been partially built up; the study of their efficacy will allow proper maintenance and/or adjustments for the upcoming years. In Greece, Italy, and Ireland OALs, this long-term efficacy study will permit the selection of the most suitable NBSs. In future climate change scenarios, these distinct classes of OALs allow the identification of the best NBS projecting and maintaining strategies.

For the sake of brevity, the detailed results of the presented models’ chains are left to further studies. Nevertheless, preliminary considerations highlight the robust link between a large portion of the presented NBSs and land-use change.

The main goals of this work are: (1) to outline the interdisciplinarity behind the study of NBSs to mitigate natural hazards; (2) to present models’ chains that allow the estimation of specific NBSs’ efficiency for each OPERANDUM OAL; (3) to show that a considerable portion of the presented NBSs is related to changes in land use.

2. Materials and Methods

The models presented in this paper are mainly well-established, and their applications cover a broad variety of science disciplines (see the OALs subsections for descriptions and references). This diversification allows us to present the models’ chains as a general approach to the NBS efficiency study. For each OAL, an intuitive flow chart describes the numerical model chain adopted: rectangles indicate the numerical models, while cylinders indicate the dataset (Section 3). Each color is associated with a type of model or dataset: the hydrologic, hydraulic, hydrodynamic, and wave models are shown with blue, the atmospheric models with light blue, the bio-geo-chemical models with green, and the solid earth models with yellow (Section 3). The novelty we propose is represented by the merging and coupling of these models with the uniform purpose of estimating the long-term efficiency of specific NBSs. We will show that this aspect related to NBSs does not require a new modeling effort. A proper association of models adopted in specific fields of interest is sufficient, also for future climate scenarios. In the following subsections, the OALs, the hazards, and the related NBSs are introduced, with a focus on the modeling activity developed to estimate NBSs’ long-term efficiency. We will present the rationale behind the selection of each model and the numerical modeling chain integrated into each OAL. In Figure 1a, the locations of all the selected European OALs are shown. Figure 1b–g is zoomed on each OAL to show the adopted models’ domains. The OALs are placed in distinct climate zones and present local peculiarities that will be described in detail in Section 3.

3. Results

3.1. OAL Finland

The OAL Finland (OAL-FI) is located in Lake Puruvesi, which represents an exceptionally oligotrophic, clearwater, and humic lake (416 km²) belonging to the Nature 2000 network (https://ec.europa.eu/environment/nature/natura2000/index_en.htm, accessed on 15 February 2021). Lake Puruvesi is surrounded mostly by boreal forests (92% of land area), which grow either on mineral or organic (peatland) soils. The main land use is forestry, with minor areas in agriculture, peat harvesting, and urban land use. Activities related to these land uses and infrequently occurring high runoff periods that might become more frequent in the future climate impose nutrient (nitrogen, phosphorus) load risk. This gradually increases the eutrophication of the lakes and threatens the recreation, fishing, and biodiversity of the area. The study area, Lake Vehkajärvi-Kuonanjärvi sub-catchment (73 km²), is located in the north-western part of the Lake Puruvesi catchment in Southern Savonia, Finland (Lat 61°58′45″ N, Lon 29°12′10″ E), Figure 1a. The two largest lakes in the study area have been classified as poor and moderate in their ecological status, and thus require intensive attention. Actions need to be taken to improve the quality of the water in the sub-catchments with poor and moderate conditions to ensure the excellent ecological condition of Lake Puruvesi in the future, following the EU water framework directive.

Nutrient and sediment loading emerges as a result of complex processes that have spatial and temporal variability. To evaluate NBSs’ effectiveness to mitigate nutrient and sediment loading, we have used (see Figure 1b for spatial domains) the National Integrated Modelling (NIM) framework, including the models RUSLE 2015 [41,42], NutSpaFHy, VEMALA [43], and LLR [44]. We have evaluated NBSs by means of two approaches: (i) including certain NBSs into the NutSpaFHy model and (ii) using possible loads simulated by the NutSpaFHy model to assess the potential of NBSs to handle load increments. In the case of OAL Finland, the NBSs considered are continuous cover forestry (CCF), riparian buffer zones, constructed wetlands, overland flow areas, sedimentation ponds and pits, peak flow control structures, and submerged dams. Riparian buffer zones are uncut areas between water bodies and tree-felling zones. By avoiding the disturbance of soil and saving the tree canopy and shrub layer, erosion and water load are reduced. Buffer zones also play an important role in sustaining biodiversity and landscape. The effects of the forest practices on the recipient watercourses can be mitigated by varying the width of the buffer zones and by following shapes of the landscape. The CCF method is a forest management regime without clearfelling. The regime involves the maintenance of a forest canopy at all times. After the felling of individual large trees, the remaining trees accelerate their growth, and new trees grow from the undergrowth reserve. Eventually, more emerge through natural regeneration. It is assumed that less nutrient and sediment leaching occur using the CCF regime compared with clear-cutting. This is because the forest is covered with vegetation all the time, and therefore the leaching of sediment and nutrients occurs less. Additionally, in peatland forests, the remaining trees keep the soil dry enough. Therefore, ditch network maintenance, which may cause additional nutrient and sediment leaching, is not needed.

The RUSLE 2015 model is used to identify areas with a high risk for erosion to cover the assessment for sediment load. NutSpaFHy is used to calculate nutrient load (nitrogen, phosphorus) from forested areas. Their outputs will be used as inputs for VEMALA, which calculates nutrient load from other land use classes and produces overall nutrient load evaluation on sub-catchment level. The LLR model is needed to evaluate the current ecological status of the lakes, and it gives the target nutrient reduction levels for which NBSs’ functioning can be compared. Data from the general circulation models MPI-ESM (representing medium equilibrium climate sensitivity) and HadGEM2 (representing high equilibrium climate sensitivity) are used. They are forced by the RCP8.5 emission scenario, downscaled by regional model RCA4, and bias-corrected [45,46]. Finally, they are used to run the models in future climate conditions (2040–2050 and 2070–2080) to evaluate the effects of increased temperature and precipitation on nutrient loads, and to assess the functioning of NBSs in such conditions.

NutSpaFHy is operated in 16 × 16 m resolution on a monthly time-step. It uses several classes of input data: forest resources, weather, topography, soil, clear-cut, and continuous cover forestry areas (see Figure 2 for the entire sequence of models. NBSs (buffers zones and continuous cover forestry) can be included by running the model with varying clear-cut and CCF scenarios. When a clear-cut occurs, the model updates the volume, height, leaf-area, and ground vegetation biomass data for the grid cells under clear-cut. Similarly, in the case of CCF, these are adjusted to examine their impact on nutrient balance. Buffer zones are considered by limiting the cut area around water elements. Outputs are reported at the sub-catchment level.

We extracted forest compartment data from the Finnish Forest Centre and the data include information about soil types, fertility, and woodstock characteristics [47]. These include age, volume, and other physical dimensions. Compartment data also include the timing and type of the next harvesting operation based on tree growth simulation and common harvesting models. The same information is available to the landowners and forest operators. We assume this data as a basis of harvesting scenarios, providing a rather realistic basis for modeled nutrient export.

The business-as-usual (BAU) harvesting scenario represents nutrient export caused by the clear-cuts, based on historical harvesting data from the area. We have processed the historical yearly average final felling area from the past ten years in the modeling area, based on forest use declarations made before harvesting to the Finnish Forest Centre. During the last ten-year period, the average total clear-cut area in the Lake Vehka-Kuonanjärvi sub-catchment area was 219 ha, of which 17% was conducted in peatland. For the BAU scenarios (Figure 3), a random selection from all potential compartments has been chosen considering the yearly averages of clear-cut areas on mineral soils and peatlands, following economic regeneration maturity. For the logging scenario including the NBSs (buffer zones and CCF), the compartment selection is conducted similarly to the BAU scenario. However, the areas identified suitable for CCF are chosen, and then to reach the same tree volume harvest clear-cut areas are also chosen. The current forest harvesting operations involve already certain buffer zones, which generally can still involve the removal of trees while not allowing soil disturbance near water bodies. However, in the logging scenario including the NBSs, the effect of wider buffer zones and avoiding tree removal can be examined. For future climate impact modeling, the periods between the years 2040–2050 and 2070–2080 are of particular interest. The uncertainties induced by making estimations on how forests will be managed during the next 100 years in the Lake Vehka-Kuonanjärvi sub-catchment area must be minimized. To this goal, we will use logging scenarios built for current conditions and run the model in 10-year periods with future climate data. This will allow us to distinguish the effect of climate change on the functioning of NBSs.

3.2. OAL Greece

The OAL Greece (OAL-GR) is located in the Spercheios River basin, in Phthiotis (central Greece, see Figure 1a). It is one of the largest rivers in Greece (80 km long) and with a drainage area of approximately 1800 km². The Spercheios river basin is mostly prone to flooding due to extreme precipitation events (flash floods) and, to a smaller extent, drought [48]. Droughts in OAL-GR are almost entirely attributed to the overuse of water sources for irrigation in specific agricultural areas (socio-economic drought). Therefore, the NBS deployed in this location is a water retention pond, to first and foremost address the issue of flooding in a specific area. The NBS is relatively small in size (especially compared to the entire length of the river), but the selected location was ideal due to several factors—the most important being the position since the overflow of the river can cause damages to infrastructures (road and train network).

Numerical solutions are ideal in the context of the OAL-GR since there are no automated hydrometric stations installed and actual streamflow data are taken sporadically. To this goal, the modeling activity for OAL-GR is focused on determining flood levels, as well as groundwater resources. To achieve these, systems at various scales have been used. The need for this multi-scale approach stems from the fact that the NBS is relatively small compared to the rest of the Spercheios River. However, it is affected by the movement of water in the entire basin.

To this end, three different models have been applied. TUFLOW [49] has been used for the flood levels and the flood mapping outputs. It consists of a suite of advanced numerical engines and supporting tools for simulating free-surface water flow for urban waterways, rivers, floodplains, estuaries, and coastlines. It can resolve water flow at very small scales, suitable for the modeling required within the above-mentioned scopes. The final aim is to provide information on the impact that the NBS has on the flood levels and the flood mapping outputs. Thus, outcomes are used to evaluate the efficiency and the advantages of the intervention. TUFLOW uses precipitation data from the ERA5 and EURO-CORDEX as input, as well as streamflow data from WRF-Hydro presented below. TUFLOW covers an area of 12 km × 10 km (Figure 1c), with a 10 m grid resolution, and it represents the main tool for examining the NBS efficiency for flood hazards. The model also uses a 10 m Digital Elevation Model (DEM) provided by the Greek Land Registry to define the river network and the surrounding topographic characteristics, and the CORINE dataset for land uses. The hydraulic model is used to calculate 10% and 20% annual exceedance probability events, before and after the NBS deployment.

WRF-Hydro [50] has been developed to facilitate an improved representation of terrestrial hydrologic processes related to the spatial redistribution of surface, subsurface, and channel waters across the land surface. Simultaneously, it aims to facilitate the coupling of hydrologic models with atmospheric models, and it is ideal for the operationalization of the NBS monitoring. In this framework, its main goal is to provide streamflow data to the hydraulic model (TUFLOW) in selected cases to better represent the incoming flow outside the bounds of the area of interest. Additionally, WRF-Hydro is more suitable for operational forecasting of river conditions and can be used for the operationalization of the NBS in the future. WRF-Hydro covers the entire Spercheios River basin (80 km × 50 km) with a 100 m resolution (Figure 1c). These data are downscaled to 10 m to be used as input in TUFLOW. WRF-Hydro uses a 90 m DEM to define the river network and global datasets for vegetation and soil use and meteorological forcing data from ERA5/EURO-CORDEX. WRF-Hydro is used in selected, high-impact cases.

Finally, the hydrological model MIKE SHE (developed by the Danish Hydraulic Institute Water and Environment—DHI Water and Environment) has been used to examine the impacts of the NBS on groundwater resources in the area. MIKE SHE is a physically based distributed model, able to simulate all hydrological domains within the land phase of the hydrological cycle in a river basin [51]. MIKE SHE is fully integrated with the channel-flow code MIKE Hydro River. This is a one-dimensional model that can simulate water flow and level, water quality, and sediment transport in rivers, flood plains, irrigation canals, reservoirs, and other inland water bodies [52]. The grid spacing for the simulations was set to 400 m. The topography of the Spercheios river basin was defined from a 50 m × 50 m DEM of the area, after resampling.

All the models can be initialized with any type of meteorological data (GFS, ERA5, EURO-Cordex). The whole models’ chain is illustrated in Figure 4. Following this procedure, we address all relevant scales for OAL-GR’s NBS resolution, starting from mesoscale and going up to 10 m.

3.3. OAL Ireland

The OAL Ireland (OAL-IE) is located in the Ringsend region in Dublin (Figure 1a), which is economically extremely important. About 10% of Ireland’s entire gross domestic product is generated from this area. The area is surrounded by River Dodder in the West and River Liffey in the North, while the Irish sea covers the eastern boundary. The Ringsend area belongs to the Dodder river basin. Since the area is surrounded by water bodies, it is prone to pluvial, fluvial, and tidal flooding.

The length of the Dodder River is approximately 27 km and the catchment area is around 113 km². Due to high slopes in the catchment area, the upper and middle section of the river is highly susceptible to flooding during periods of extreme rainfall events. Recent flood events prove that the river can exhibit a maximum flow up to 250 m³/s. As the river exhibits a small basin with steep slopes, it usually takes 2–3 h to reach the water from upstream at the mountains to downstream at Ringsend, potentially creating fluvial flooding (flash flooding). Pluvial floods are common as a consequence of extreme rainfall in such a lower elevation area. Controversially, tidal flooding was found to be rare in the OAL-IE surroundings. Hence, this part of the study focuses on the investigation of the changes in both pluvial and fluvial flooding in OAL-IE due to climate change, and the identification of the effectiveness of a hypothetical nature-based solution (NBS) to reduce floods in the area.

The Soil Water Assessment Tool (SWAT) model has been used to simulate streamflow at River Dodder for understanding fluvial flooding and pluvial flooding in terms of rainfall. The advantage of the SWAT model is that it accounts for the fluvial as well as the pluvial flooding in the area. It is a physically based, semi-distributed rainfall-runoff model. It requires as inputs land cover, soil, and five meteorological variables’ (rainfall, daily maximum and minimum temperature, relative humidity, solar radiation, and wind speed) information. Several pieces of research have been conducted using the SWAT model. A detailed set of literature on SWAT can be obtained via https://www.card.iastate.edu/swat_articles/ (accessed on 14 April 2021). The incorporation of NBS implementation in the study area can be simulated by changing the local land cover and soil properties. The effects of future climate change can be estimated by feeding the future projected values of the above-mentioned meteorological variables.

These future projections can be obtained from a set of Regional Climate Models (RCMs), corresponding to different Representative Concentration Pathway (RCP) scenarios. In this study, several RCMs were considered and two RCPs (4.5 and 8.5) representing different climatic conditions were selected to obtain future projections of the meteorological variables. The RCM outputs were generated as a gridded scale averaged over an area. In this study, around 12.25 km spacing was considered for every grid. The RCM covered the entire Dodder river basin (Figure 1d). Bias correction of the projected future meteorological variables was performed using the Quantile regression-based bias correction approach.

The SWAT model was developed for the entire Dodder river basin (Figure 1d), as the Ringsend is located at the junction of River Dodder and River Liffey. The basin boundaries were obtained from the Environmental Protection Agency, Ireland website. To develop the physical-based model, DEMs were obtained from COPERNICUS EU-DEM (25-m resolution) and Bluesky DEM (1-m resolution). The land cover and soil maps of the watershed needed to be included in the SWAT model once the watershed and the river network were formed. Those maps are essential for the estimation of the water movement over the land surface, the estimation of evapotranspiration from the surface, the amount of water percolation through seepage, the amount of subsurface flow, and the groundwater recharge. In this study, the land cover map was obtained from the CORINE land use map, available for Europe (https://land.copernicus.eu/pan-european/corine-land-cover/clc2018, accessed on 14 April 2021). The land cover classes do not match with the required land cover classes for the SWAT model. Thus, the land cover classes from CORINE data were reclassified to the following seven classes that are present for the Dodder river basin: Agricultural land, Barren land, Forest area, Pasture, Urban area, Waterbodies, and Wetlands. The soil maps for Ireland were obtained from the General Soil Map of Ireland, available here http://gis.teagasc.ie/soils/downloads.php (accessed on 14 April 2021). Similar to land cover classes, the soil classes also needed to be reclassified for the SWAT model. The slope is classified into five classes: 0–2%, 2–5%, 5–10%, 10–25% and >25%, where an increase in the value means steeper slopes. Based on all the information, the catchment was divided into multiple sub-basins, which were then further subdivided into hydrological response units (HRUs) using SWAT.

The observed values of the needed meteorological variables were obtained from Met Eireann and Dublin City Council’s monitoring stations as well as from newly installed weather stations through the OPERANDUM project. The HRUs were the smallest scale sub-basins that consist of homogeneous land-use, management (slope), and soil characteristics. Once the HRUs for the watershed were formed, SWAT simulates the flow dynamics in the watershed. The water balance of each HRU is represented by four storage volumes including snow, soil profile (0–2 m), shallow aquifer (typically 2–20 m), and deep aquifer (>20 m).

In the SWAT model, lateral subsurface flow in the soil profile (0–2 m) is calculated simultaneously with percolation. A kinematic storage routing is used in SWAT to predict the lateral flow in each soil layer. The procedure is based on the degree of slope, slope length, and saturated hydraulic conductivity of the soil layer. Lateral flow occurs when the storage in any layer exceeds field capacity after percolation. Groundwater flow contribution to total streamflow is simulated by creating shallow aquifer storage. Percolation from the bottom of the root zone is considered as a recharge to the shallow aquifer. In order to simulate the surface runoff, the water balance equation requires information of five meteorological variables at daily scale: rainfall, mean temperature (which is considered to be the average of the maximum and minimum temperature in the study), wind speed, relative humidity, and solar radiation. To reduce runoff that contributes to flooding, the conversion of car parking to natural grassland was considered to be the hypothetical (potential) NBS in the Ringsend region. The general flood analysis framework by combining different models (RCM, SWAT) and data (basin characteristics, meteorological data) is shown in Figure 5.

3.4. OAL Italy

The OAL Italy (OAL-IT, see Figure 1a) is composed of different sites (Bellocchio, Panaro, and Po di Goro), involving distinct categories of natural hazards. The Bellocchio site is located in the coastal area of the Emilia-Romagna region at Bellocchio lagoon (North of Italy, see Figure 1e). The area is subject to coastal inundation and erosion due to storm surge and waves. Two kinds of NBSs proposed to mitigate these hazards are sand dunes and marine seagrass set up in the coast’s proximity. Concerning artificial dunes, reference [53] documented that along the high-energy Oregon coast, dunes plus cobble grounds have provided an acceptable level of protection to the areas behind the NBS. As regards the marine seagrass, several authors documented the potential to decrease the energy of the waves and the intensity of the storm surge (e.g., [54]). The combination of the above-mentioned NBSs (namely NBS1 and NBS4 for sand dunes and marine seagrass, respectively) could represent a promising and innovative approach to coastal protection solutions.

Both NBSs will be assessed through a three phases process illustrated in Figure 6. The NBS1 is implemented following the present legal constraints, while numerical models are implemented in the region of interest for NBS4. Secondly, the numerical simulations are carried out in the present and future climate scenarios with and without NBSs. Eventually, an analysis of the model simulations is carried out to assess the fitness for purpose of the two NBSs.

In OAL-IT referring to Bellocchio, the estimation of the efficiency of the NBSs is intricate due to the considerable space-time resolution required. In addition, single numerical models are not able to incorporate all the spectrum of the processes involved. Thus, a multi-model strategy must be preferred. Storm surge can be modeled by means of two numerical models: (1) a sea-level model (considering astronomical tides, atmospheric pressure, and winds), and (2) a wave model (simulating wind waves). Eventually, a third model for coastal erosion considers sediment transport, morphological evolution, and bathymetry changes.

The sea-level model is based on the hydrostatic primitive equations for the ocean and accounts for thermodynamic and hydrodynamic processes. The horizontal and vertical scales are in the order of hundreds of meters and tens of meters, respectively. In particular, the SHYFEM unstructured grid hydrodynamical model [55,56,57] is used here to evaluate the fitness for purpose of submerged aquatic vegetation on the sea level extremes, under present and future climate conditions.

The unstructured WAVEWATCH III (WW3) model has been chosen to simulate the wave characteristics in the Emilia-Romagna coastal strip. WW3 is a community wave modeling framework that includes the latest scientific advancements [58] in the field of wind-wave modeling and dynamics. The wind-wave model integrates a balance equation for the wave action density spectrum. This model solves for the significant wave height, the wave direction, the wave mean and peak period and the Stokes velocities induced by the wave motion. It is designed for resolutions in the order of a few hundred meters, as for the hydrodynamic current model. Then, the Stokes velocities can be coupled to the hydrodynamic model to produce more realistic currents at the surface. These are used in advanced formulations of bottom stress given by the combination of current and bottom orbital velocities.

Eventually, the third model estimates the coastal erosion, given the sea level and the waves coming from the first two models at large scales. The spatial scale of the model is in the order of k-meters. The code contains a depth-averaged advection-diffusion sediment transport equation with a source-sink term, based on equilibrium sediment concentrations. The model is therefore able to simulate the effects of storm surge on the dunes by including the avalanching process. In order to assess the NBS reliability of OAL-IT Bellocchio, consisting of a naturally reinforced dune, the morphological model XBeach has been used. The model allows the simulation of many and complex near-shore hydrodynamic and morphodynamic processes at local scales. It has been applied and validated for storm impacts on the dune and urbanized coasts for dune safety assessments. Thus, it is suitable for this type of application. For the hydrodynamic processes, the model includes short wave transformation (refraction, shoaling, and breaking), long-wave transformation (generation, propagation, and dissipation), and flow [59]. Originally developed to simulate impacts on sandy coasts with a domain size of kilometers and on the time scale of storms, the model has then been applied to other types of purposes and for long-term prediction [60,61]. A correct representation of the bathymetry is crucial for modeling morphological processes. In this study, the bathymetric elevations were derived from the merging of topographic and bathymetric surveys. A local topographic survey performed specifically for the OPERANDUM Project has been used to generate a detailed topography of the natural dune and the area where the NBSs are implemented.

The modeling strategy is schematically illustrated in Figure 7. The three models will be coupled offline for two scenarios: the first one, related to the present climate, and the second one, to the future climate considering RCP 8.5 (emissions in RCP 8.5 continue to grow until 2100). The ocean and atmospheric forcing for the present climate are derived from ECMWF analyses and CMEMS data (https://marine.copernicus.eu/, accessed on 14 April 2021). For the present and future climate, the ocean and atmospheric forcing will be extracted from the MedCORDEX data set [62]. The ECMWF-CMEMS driven simulations for the present time can set the accuracy of the model simulations for the sea-level, wave, and coastal erosion conditions in comparison with in situ data.

The OAL-IT Panaro river site (Figure 1f) is devoted to the modeling of the dynamical interaction between water, soil, and the NBSs during extreme flooding events. In particular, the NBS conceived to mitigate the risk of riverine inundations refers to the establishment of a dense vegetation cover on the earth embankment intending to reduce the likelihood of its failure due to erosion phenomenon. This target could be achieved using deep rooting perennial herbaceous species, whose aerial parts are expected to reduce the mechanical impact of raindrops, while belowground parts (deep thin roots) are demanded to mechanically reinforce the embankment-forming materials, facilitate drainage in the topmost layers of the earth embankment, and promote the plant water uptake. OAL-IT Panaro relies on experimental planting and data collection activity carried out along a portion of the internal slope of the river embankment. However, testing the NBS for present and future hydraulic loads, as well as over a considerable space-time resolution (i.e., extensive embankment system and plant growing period), is challenging, thus requiring a complex modeling framework based on multiple models.

The modeling cascade ideated for this scope is shown in Figure 8. The analysis is performed by means of rainfall-runoff and hydraulic modeling at a larger (catchment and river) scale and employing advanced CFD (Computational Fluid Dynamics) approaches at smaller scales (river stretch). The TUW rainfall-runoff model (a semi-distributed hydrological model [63]) can generate the continuous hourly time-series of streamflow expected at the closure section of the river catchment, corresponding to the retention basin location, in case of different meteorological input. The analysis will consider both historical events, based on the observations recorded in a set of temperature and rain gages managed by the regional meteo-hydrological service (SIMC-ARPAE), and future climate scenarios. The bias correction of climate scenarios at a daily scale for precipitation and temperature will be carried out using the scaled distribution mapping (SDM) procedure [64]; precipitation adjusted values will then be post-processed to retrieve assessments at a sub-daily scale by using the Random Parameter Bartlett–Lewis (RPBL) model as suggested in [65]. The propagation of the floods estimated from the TUW model, from the closure section of the catchment up to the section considered for the OAL-IT Panaro, is performed employing the HEC-RAS hydraulic model. Based on the UNET code for solving the De Saint–Venant equations [66], HEC-RAS is used to simulate hydraulic loads at the OAL site in case of design flood events (events at given return periods) or long series of flowing conditions associated with potential climate change scenarios. Simulation outputs are expressed in terms of hydraulic variables, such as water depth and mean velocity, flood duration, velocity of flood rising and lowering phases. Finally, outcomes of the numerical modeling will serve as boundary conditions for further additional investigations. CFD models, developed within OpenFOAM environments, are used to estimate the interaction of water and NBS during the flood and the fluid–structure interaction. This is to obtain the forces acting on the NBS and the modification of the NBS shape caused by the action of the fluid on the embankment.

Different cases of river flooding are compared, with flow conditions obtained by upper scales models (TUW + HEC-RAS) and used as boundary conditions for the CFD analysis. The geometry of the river is obtained by an experimental campaign performed on the site. The NBS is modeled as a solid that can be modified in shape as a result of the fluid effect (employing fluid–structure interactions modeling). These results allow the characterization of the NBS capacity to withstand an extreme event such as flooding.

This OAL-IT Po di Goro site (Figure 1e) is devoted to the modeling of the salt wedge intrusion along a branch of the Po river delta. The estuary of the Po di Goro branch is a well-known “salt-wedge” estuary: the interannual average of the salt wedge intrusion is around 15 km (according to ArpaE monitoring campaigns over 2003–2017), but the entering salt ocean water may advance more than 20 km upstream during low runoff. The modeling strategy is schematically illustrated in Figure 9. The core of the system is a 2-layer 1D Estuary Box Model, the so-called CMCC EBM [67]. It aims to (i) represent the net river release at the estuary mouth in terms of volume flux and salt flux and (ii) estimate the length of the salt-wedge intrusion. The river runoff at the estuary head and the ocean inflow (baroclinic in flow through the bottom and barotropic inflow during the flood tide phase) at the estuary mouth are the inputs needed by the CMCC EBM model. In the framework of the Operandum project, the CMCC EBM has been applied to the Po di Goro estuary to produce both historical simulations (1982–2011) and middle-term scenario (2021–2050 under RCP 8.5) of both the salt wedge length and the salinity at the river mouth. The CMCC EBM preliminary results clearly show a stronger intrusion of the salt wedge in the middle term. The historical simulations and the scenarios performed with the CMCC EBM for the Po di Goro estuary indicate that the foreseen increase in the salt wedge length in the middle-term is about 14% (up to 50% in summer). Moreover, the projected increase in the salinity at the river mouth is about 9% (up to 35% in summer). The average, the minimum, and the maximum values of salinity at the Po di Goro mouth projected by the CMCC EBM are considered as reference values to investigate the behavior of two halophyte species selected as the site-specific NBS. These plants are expected to effectively reduce the saline intrusion owing to their high salinity absorption capacity.

3.5. OAL UK

The OAL United-Kingdom (OAL-UK) is located adjacent to Catterline Bay in Northeast Scotland (Aberdeenshire; WGS84 Long: −2.21 Lat: 56.90; Figure 1a). The site belongs in the climatic temperate humid region, with a mean annual temperature of 8 °C and mean annual rainfall of ca. 600 mm [68]. The OAL’s landscape is dominated by sloped terrain and cliffs running into a cobbled beach, which encloses the sea in a half-circle bay where seals and seabirds breed and shelter. The slopes and cliffs are mainly vegetated with early successional plant communities (i.e., herbs and grasses) surrounded by a mixture of areas dominated by riparian trees and shrubs. Severe shallow landslides, surface erosion, and coastal erosion episodes have occurred in the past because of prolonged rainfall periods and sea storms coincident with spring tides.

Multiple plant-based NBSs are planned at OAL-UK to mitigate and manage shallow landslides and erosion. Examples of these NBSs are live cribwalls and high-density plantations of woody vegetation. A live cribwall is a gravity wall made of a cellular structure with logs and living stakes or plants in containers, with the objective that the future development of the plant substitutes the log structure following its natural degradation or decay. This intervention is used as a retention wall in the stabilization of slopes with gradients up to 60°. The cribwall structure is made of softwood timber logs (e.g., logs of peeled conifer) deployed perpendicularly and horizontally with respect to the slope. The space between the logs is filled with compacted earth materials and living stakes, and cuttings or plants in containers are planted on the front part in high density to retain and reinforce the soil materials. Plant fascines can also be deployed within the front part, which will contribute to retain moisture and foster the establishment of vegetation within the cribwall. The structural logs are secured with nails or steel bars to ensure that the cribwall structure works as a block.

The high-density plantation of woody vegetation is a simple NBS technique that can substantially contribute to the protection and reinforcement of slopes by using plant cuttings and stakes from local woody vegetation with the ability to propagate from stems. To provide effective slope reinforcement, it is important to establish a high-density plantation at the toe of the slope. However, planting location can be selected according to environmental features within a given slope, e.g., dryness, aspect. Woody vegetation will provide mechanical and hydrological stability to the slope forming materials and it will provide a ground cover protecting the slope against surface erosion.

The modeling activities at OAL-UK aim to evaluate the performance of plant-based NBSs against shallow landslides and erosion under current and future climate change scenarios. The modeling activities are distributed into two nested spatial scales, i.e., slope and OAL scales, for which several model systems are used. All the models are evaluated in a daily time-step and use forcing meteorological data for the current and future climate conditions originated from the HadGEM2-ES model for the RCP 8.5 scenario.

At the slope scale, the combined reinforcement performance of cribwall structures and live vegetation is assessed with a novel numerical model operating under changing ground cover and meteorological conditions at a daily time step—i.e., the LiveCW model. LiveCW ultimately evaluates slope stability under changing soil moisture conditions by calculating the Factor of Safety (FoS) with a limit equilibrium method (LEM), which employs Bishop’s ordinary method of slices [69]. The geometry of the slope and potential slip surfaces are established based on 2D topographical profiles derived from a 2 m resolution digital terrain model (DTM) using basic circumference and linear geometry. The combination of basic approaches on circumference geometry, geotechnical design, and structural mechanics is employed to estimate the cribwall geo-grid and face-wall contribution to slope stability through internal compound stability (ICS) analysis. The service life of the cribwall structure is considered by incorporating the effect of timber decay on soil reinforcement under changing climatic conditions using the model proposed by [70]. The hydro-mechanical reinforcement effect of vegetation on the cribwall and the engineered slope is assessed under changing meteorological conditions following the establishment of the slope and cribwall geometries. To this end, the workflow established in the Plant-Best model [71] is incorporated into a novel plant-slice module (PSM). Plant-Best uses plant, soil, and meteorological inputs to estimate the mechanical and hydrological effect of vegetation on soil reinforcement at interoperable spatial scales. PSM splits the simulated slope profile into multiple slices and estimates the plant-soil, hydro-mechanical reinforcement [72] within each slice before coupling this effect with ICS. LiveCW is written in the open-source, statistical language R, making the model versatile and customizable. ICS and LEM approaches are widely used in geotechnical design (e.g., [73]), yet models evaluating the stability effect of timber-based cribwalls are scarce and only established at the design level [74]. To the best of our knowledge, the combination of the reinforcement effect of cribwalls and vegetation into one model has not been attempted yet. Additionally, it sets the basis for building numerical models assessing other types of timber and plant-based NBSs deployed in OAL-UK such as live lattices or gratings. Although LiveCW is built in a 2D fashion and at the slope scale, the incorporation of surface-based variables from the plant, soil, and atmospheric compartments will facilitate its upscale and reproduction into alternative environmental settings. The whole model strategy for OAL-UK is shown in Figure 10.

4. Discussion and Conclusions

The natural hazards mitigation policies are gaining a leading role in the climate change adaptation framework and disaster risk reduction [75,76,77]. Herein, we have presented a series of specific sites where distinct hydro-meteorological hazards are coped with NBSs. The means of OALs is therefore gaining a wider significance. They allow the cooperation of a variety of local stakeholders and the scientific community to find the best solution for the mitigation of specific natural hazards. There is a strong need for this kind of cooperation to allow the decision-makers to have a broad significance of the action to be undertaken.

Eventually, OALs serve to demonstrate that the modeling strategies to estimate NBSs’ long-term efficiency can follow a similar path, despite the local scope. The multidisciplinary approach the NBSs are required for is glaring in the models’ chains presented.

In all the given cases, the climate models’ outputs serve as inputs for all the future scenarios modeled. This general procedure allows the identification of the long-term key-role physical parameters in the estimation of NBSs’ efficiency. This is possible through the analysis of the current (where available) status of the NBSs. If we consider the example of the cribwall at OAL-UK, the monitoring variable is the factor of safety. Thus, the climate projections of interest are related to the variables controlling the FoS (e.g., groundwater alterations related to changes in precipitation trends). At this stage, the major challenge in the NBSs’ case is represented by the needed space scale. NBSs are generally built up to mitigate local hazards, in the range of hundreds or thousands of meters. On the other hand, climate models’ scales are in the order of tens of kilometers. Thus, the implications of future climate scenarios in the efficiency of NBSs can be considered only by systematic downscaling of the variables of interest. Downscaling techniques have been widely studied for a multitude of environments and purposes (e.g., [78,79,80,81]). Nevertheless, it is hard to detect common strategies, especially in locations where local effects (e.g., strong vertical gradients, limited lakes) can strongly influence the climate. Hence, the choice of the most appropriate downscaling technique depends on the available information in the region where the NBS is located. This aspect is not the target of this study, but the models’ chains presented here emphasize the need for a systematic study of the region before projecting NBSs. This aspect represents the main limitation of the presented models’ chains since the climate input data have already been biased. Additionally, the downscaling techniques have a strong dependence on the available sets of data in the area of interest. Anyhow, this limitation is implicitly correlated to the current inadequacy of climate models to catch lower-scale phenomena.

The application of the presented models’ chain to the above-mentioned NBSs is under development. Nevertheless, the available preliminary results demonstrate that models are able to catch the NBSs’ behavior in future climate scenarios.

NBSs will represent an inclusive tool related to land management. For this reason, we believe that the administration of natural resources will also be influenced by the presence of NBSs. If we consider the OAL-FI and OAL-IT, it is clear how the environments show considerable differences. While in OAL-FI continuous cover forestry is adopted to control the nutrient loading, in the OAL-IT Po di Goro site the specific vegetation will control the salt intrusion. Despite the distinct purpose of the NBSs and related modeling activity, a common finding is glaring: the adoption of land management to cope with natural hazards. Indeed, other presented NBSs show a more immediate effect. If we look at the OAL-UK, live cribwall installations prove to be effective in the control of landslides. However, in this case, the vegetation can have a crucial effect on the long-term stabilization of the landslide. Additionally, in the OAL-IE site, a land usage change has been proposed as an NBS. Eventually, in the OAL-GR, a water retention pond has been chosen as an NBS, proving again the importance of land management. Unsurprisingly, this aspect is not always evident in the models shown here since they involve multiple sections of the OALs simultaneously. Indeed, land changes can be studied formally (e.g., [82]) or by employing cost–benefit analysis (e.g., [83]). In both cases, the long-term effects of the processes are generally the focus of the studies. In this view, the formal approach developed in the OPERANDUM project for the estimation of NBSs’ long-term efficiency shares crucial aspects with the land management policies. The implications of these results could have a relevant impact on the natural hazards mitigation policies. This is also evident in recent works about the monitoring and modeling of NBSs [84,85]. OAL and NBS concepts could provide future standards for these policies, delivering easy-to-handle information to decision-makers and stakeholders. Therefore, we recommend the implementation of NBSs in other locations to also face different hazards. This will allow comparisons with the presented approach. Extensive results in numerous locations will allow a systematic generalization of the OAL-NBS implementation.

The models’ chains presented herein represent an important step in the conceptualization of NBSs in an environmental framework. Future developments will further show how specific kinds of NBSs can have direct consequences in land and natural resources management.