Combining machine learning and spatial data processing techniques for allocation of large-scale nature-based solutions

Hydro-meteorological hazards such as floods, landslides, droughts, and heatwaves have increased globally due to human activities and climate change (Singh et al. 2016). These hazards produce environmental damage, life losses, and world economy stress every year (IPCC 2014). To reduce impacts of climate change, nature-based solutions (NBSs) are suitable measures, especially considering the conventional hydraulic structures are not sustainable and adaptive. NBSs are solutions and actions that mimic nature, integrating hydrology, geomorphology, hydraulic, and ecological dynamics (European Commission Directorate-General for Research & Innovation n.d.; Haring et al. 2021). NBSs can offer multiple benefits through multi functionality, but their effectiveness varies due to the scale of applications.

NBSs can be categorised into three types based on the scales of analysis: fine, local, and regional scales (Somarakis et al. 2019). Fine-scale NBSs include gardens, small parks and vegetated roofs, and walls. They are usually implemented to mitigate noise and heat islands, support biodiversity, reduce flood risk, and save energy. Local-scale NBSs are larger parks, forests, and other green spaces that impact the local environment. These measures manage the local environment, enhance air quality, protect wildlife, and manage water. Regional-scale NBSs are applied in rural, highland and seaside regions, such as swamps and river restoration, dunes, and wetlands. Fine and local-scale NBSs are also called small-scale NBSs, while regional-scale measures are considered as large-scale NBSs. Large-scale NBSs provide more co-benefits than the small-scale ones, considering they are applied in larger areas and interacting more with catchment processes (Watkin et al. 2019; Ruangpan et al. 2020; Kumar et al. 2021; Vojinovic et al. 2021). The co-benefits include long-term solutions, job generation, carbon sequestration, and enhancing ecosystem and biodiversity.

Examples of large-scale NBSs include rainwater harvesting (WH), wetland restoration (WR), and natural riverbank stabilisation (NRS). WH involves intentionally collecting rainwater from surfaces such as roofs, ground surfaces, and watercourses, and storing it in structures or the soil profile (Malesu et al. 2006). This practice effectively mitigates the impacts of floods and droughts in rural areas. WR focuses on identifying degraded areas within wetlands and restoring them to their original state, thus rehabilitating hydrological functions (World Bank 2008). Wetland degradation often results from land use changes to croplands. Wetlands play a crucial role in reducing flood risks by attenuating flow peaks and mitigating droughts through water storage. NRS are measures implemented to manage erosion and protect and stabilise riverbanks (Allan & Castillo 2007). Rivers undergo constant geomorphic processes, including erosion, sedimentation, and bank instability, particularly during extreme events. NRS implementation involves bioengineering and green solutions such as riverbank revegetation, reducing the impacts of floods and landslides.

Implementation of NBS projects involves key phases such as scoping, planning, design and construction and operation and maintenance, and in some cases, decommissioning. Planning includes hazard data collection, NBS assessment and socioeconomic analysis. Locating suitable NBS sites is a crucial step at the planning phase. A systematic review by Mubeen et al. (2021) identified seven tools/methods for spatial NBS allocation. Most of the tools/methods focus on small-scale urban NBSs, except for some scholars (Guerrero et al. 2018; Mubeen et al. 2021) who developed a methodology using ‘floodplain hydro morphological landscape units’ to identify potential areas for large-scale NBSs. Following the review, Mubeen et al. (2021) developed suitability mapping methodologies for detention basins, retention ponds, floodplain restoration, river widening, forest buffers and afforestation. They combined criteria based on slope, soil type, distance from streams, land use type and road buffers to create general suitability maps, which were further refined with NBS-specific criteria. These methodologies were implemented using Esri ArcMap software's model builder functionality. However, there is still the need to develop mapping methodologies for other large-scale NBSs or to further improve the existing ones.

This article further contributes to NBS planning by developing suitability mapping methodologies for WH, WR, and NRS. For WH, Mubeen et al. (2021) developed a methodology for allocating retention ponds, which requires consideration of climatic conditions for rainwater availability. Droogers & van Kempen (2014) proposed a parametric approach that uses seven spatial datasets of which some are specific to certain WH designs. The present research does not include those datasets while considering other relevant, spatial information used as constraints. As for WR and NRS, there is no standard procedure or study for mapping suitable locations. Hence, our work provides contribution in this direction. The methods developed are tested and analysed on at least one study area where each of the three large-scale NBSs were evaluated. The study was conducted within the EC-funded H2020 RECONECT project (http://www.reconect.eu/).

As the three NBSs have different characteristics and can be applied in different settings, the datasets needed and methodologies developed are also different. The upcoming sections describe the three methodologies developed. In all cases, the ModelBuider functionality in ArcGIS has been used to create geoprocessing models (toolboxes). ModelBuilder (ESRI 2022) is a visual programming language that allows automating and documenting spatial analysis and data management processes. The tool enables data integration, analysis, and visualisation, facilitating workflow sequencing and model sharing.

Part of the common input data for the toolboxes concerning WR and NRS requires a land cover raster dataset in order to produce results with higher accuracies. Several land cover products are available, such as European Spatial Agency (ESA) WorldCover (2020), and ESRI (2020) Global Land Cover Map (ESRI), which offer a spatial resolution of 10 m. Alternatively, Copernicus Corine Land Cover (CLC), and Copernicus Global Land Cover (GLC) provide a spatial resolution of 100 m. However, it has been observed that some of these products exhibit limitations in accurately identifying certain land cover types, particularly differentiating between croplands and grasslands. This issue arises due to the large scale at which these products are developed, often covering vast areas worldwide. On the other hand, while some products can correctly identify croplands, they may contain incomplete sections due to their coarse spatial resolution. Table 1 presents a summary of the key features of the four existing land cover products mentioned above, all of which are freely accessible online.

To obtain a land cover dataset with higher accuracy specifically tailored for the study areas, one of the machine learning (ML) techniques, a Random Forest (RF) classifier was employed in a supervised classification workflow implemented in Google Earth Engine (GEE), a cloud-based platform. The RF classifier was chosen due to its ability to handle complex classification problems and provide reliable results (Rodriguez-Galiano et al. 2012; Collins et al. 2020).

The classification process primarily utilised the Sentinel Multispectral image collection, which offers a spatial resolution of 10 m. In some cases, Landsat imagery was used for time periods predating the availability of Sentinel images. The selection of images was based on their availability within the desired time range.

By computing the Kappa index, the accuracy of the RF classifier can be quantified, providing an objective measure of its performance in classifying the land cover in the study areas. Landis & Koch (1977) provide the categorised interpretation of the Kappa coefficient, which indicates the level of agreement between the classification results and the reference data. Kappa values greater than 80% and less than 40% indicate strong and low agreement, respectively, between the classification and the reference data, while values between 40 and 80% represent a moderate level of agreement (Congalton & Green 2008). This classification scheme enables the assessment of the accuracy of the RF classifier based on the calculated Kappa index.

The methodology is based on an open pond WH scheme. Considering the use of WH for diverting and collecting stormwater to prevent flooding and using stored water in dry periods, the present research uses six main input datasets – total annual precipitation, coefficient of variation of the seasonal monthly precipitation, aridity index, runoff coefficient, slope, and population density. According to Droogers & van Kempen (2014), the total annual precipitation provides vital information on the amount of water available in a basin throughout the year. The variation in the monthly precipitation is used to identify the peaks and extreme events during the rainy season. The aridity index shows which areas need WH most while the runoff coefficient indicates the areas that generate higher surface runoff. Slope, which is extracted from a digital elevation model (DEM), shows the influence of topography on the surface runoff flow. As WH is implemented to benefit communities, population density identifies the spots with higher concentration of settlements in a basin.

The input data are processed in raster format and have the same spatial resolution and projection. However, as the data have different units of measure, they must be normalised, from 0 to 100. In this process the highest value of the raster is represented by 100 and the lowest by 0. Once map algebra is applied, the result is a dimensionless raster where the highest values are the most suitable to locate WR than the lowest values.

Not all areas that satisfy the above input conditions can be used to implement WH schemes. The places where WH cannot be applied should be removed at the end. The constraints considered in this study are: transport lines (in this case, roads), water bodies (in this case, rivers) and buildings. Roads and rivers are usually available in vectors and should be converted to raster, whereas buildings are extracted from land use data in raster format. The raster cells that contain buildings, buffered roads and rivers should be assigned null values to avoid using them as potential suitable areas.

Based on meteorological and remote sensing data, the present research applies various methods and datasets to analyse and generate data. The coefficient of variation of seasonal monthly precipitation was calculated based on data statistics from rain gauge stations within or near the study area. This analysis considered the seasonal variation during the wet period. The data were then interpolated in ArcGIS to produce a raster. Annual precipitation, runoff coefficient and aridity data can be generated for any region worldwide using information from meteorological services or satellite data sources such as CHIRPS, TRMM, and PERSIANN. GEE platform was employed to obtain raster data for annual precipitation, monthly seasonal variation coefficient (Banerjee et al. 2020), runoff coefficient and aridity index maps.

A DEM was used to delineate the catchment area of the study region and determine the slope. In this study, DEM Alos Palsar (2023) was used due to its higher spatial resolution (12.5 m) and the ability to delineate a basin and river network with greater precision. It is important to note that while this DEM covers most regions in South America and is freely accessible online, it is not available worldwide. For other parts of the world, other DEMs with different resolutions can be used. Population density data are usually accessed from a country's statistics office.

The facilities/buildings layer was derived from the land cover data provided by ESA. River and road networks from OpenStreetMap were buffered and converted to raster format. These buffers represent constraint areas where WH cannot be implemented and are assigned null pixel values. The buffer distance for roads follows the standard set by the American Association of State Highway and Transportation Officials (AASHTO 2020), which designates a 50-m width on each side of the road centreline as the ‘right-of-way’ area. Similarly, the same criteria of a 50-m buffer on each side of the river or water body centreline were adopted.

The international convention Ramsar, established in 1971 with the mission of preserving and promoting the sustainable global development of wetlands, provides online free access to a comprehensive database on wetlands worldwide (RAMSAR 2010). This database serves as a valuable resource for the study area, and the following methodology involves identifying areas within Ramsar wetlands where natural land cover has been converted into croplands. Both vector and raster data have been used for this purpose.

The required vector data for this analysis are the Ramsar Polygons (RAMSAR 2010), which define the boundaries of the wetland. In addition, vector data from OpenStreetMap were employed to identify railways and roads, which serve as constraints in the analysis. The buildings layer used in the analysis was derived from the ESA WorldCover (2020) product, chosen for its high accuracy in recognising buildings. This product incorporates samples from various locations worldwide, enhancing its ability to identify buildings accurately. The land cover is the result of the RF classification computed in GEE, discussed in Section 2.1. The RF classifier's accuracy relies on the sampling process, which may be limited in small, specific areas.

The toolbox workflow begins by clipping the land cover dataset, buildings layer roads and railways using the wetland polygon as a mask. This retains only the relevant portions within the wetland area for analysis. The land cover classification output is then processed to extract water bodies and crop layers, which are converted into vector format for easier integration and analysis. A 50-m buffer is applied around road and railway features based on standard guidelines (AASHTO 2020), representing the designated right-of-way area. These layers are then intersected with the crop and building layers to identify intersections that may impact WR suitability.

To refine the analysis, the resulting intersection layer is further intersected with the study area polygon, excluding the regions outside the wetland study boundary. This step ensures that only areas within the wetland study area are considered for WR measures, excluding regions beyond the study boundary. The toolbox generates a vector file indicating suitable areas for WR.

Some of the datasets required to map suitable sites for NRS are land cover information derived from the RF classification, slope data from a DEM, soil type, and erosion index. A hydrodynamic model was developed to generate the maximum depth and velocity of a river. Velocity thresholds are defined to distinguish between permissible and potentially dangerous velocities, which are dependent on the riverbank (Bureau of Reclamation 2015). For example, the threshold of permissible velocity for sandy loam soil is 0.53 m/s.

The suitability mapping process began with resampling and clipping the erosion index raster. At the same time, land cover classification rasters were used to extract four Boolean rasters (contain values true or false). Layer classified as crops and global erosion index were intersected, identifying the weakest zones. Values higher than a permissible velocity of 0.53 m/s are considered erosive velocities in the riverbanks. This value was used to create a new maximum velocity raster, which was intersected with the previous erosion index-crops raster. This is the first intermediate raster output. The normalised classification raster was combined to remove all the features not located close to the river borders (layer of sediments, water-riverbed). Separately, the slope data was intersected with the maximum depth raster, a second intermediate raster output. Finally, the two intermediate raster outputs were intersected to provide areas prone to destabilisation and require NRS. NRS are typically applied along a river, near settlements and croplands, to protect them from floods. Implementing these measures along the entire river would disrupt natural processes and lead to undesirable impacts and ecosystem disservices (Somarakis et al. 2019).

In this study, we applied the suitability mapping methodologies to four areas. However, the methodologies can be applied to any case study given data availability, as the toolboxes are generic. Below, we described the four study areas.

Datasets required to develop the toolbox were obtained from the latest water balance report by MMAYA et al. (2018), including total annual precipitation, total annual evapotranspiration, and runoff coefficient. The aridity index is a sub-product calculated by dividing total annual precipitation by total annual potential evapotranspiration. All rasters have a spatial resolution of 30 m except for precipitation, which has a coarse spatial resolution of 5,000 m. To ensure uniformity, the spatial resolutions of all input rasters and sub-products were set to 30 m.

Precipitation data from rain gauge stations in the study area were used to calculate the coefficient of variation of the seasonal monthly precipitation. The wet season was analysed separately to identify sudden heavy rain events. The coefficient of variation values were interpolated to create a raster dataset. Population data from the GeoBolivia website were interpolated and converted into raster format. Alos Palsar DEM was another source of information applied to obtain sub-products such as slope, hillside, and the delineation of the study basin. The last raster source of information was the land cover ESA 2020/ESRI products, which was used to extract the buildings layer and identify constraining locations.

To test the WR toolbox, which identifies areas with a land use change inside a wetland, two RAMSAR polygons within the RECONECT project were chosen. Figure 1(b) shows the first study area, a Ramsar wetland called National Park Biesbosch, located in Dordrecht, the Netherlands. The area of the wetland is about 103 km². The input data and general methodology depend on the availability of information and the proper characteristics of the place. Figure 1(c) displays the second study area, a Ramsar wetland called Peštersko Polje, located in Serbia. The wetland has an area of about 34 km², and it belongs to the Drina River basin.

The study area selected to test the NRS toolbox is located along the San Juan del Oro River, near Arteza town in Tarija, Bolivia. The region primarily relies on agriculture, and the study area is susceptible to unstable riverbanks and flooding, impacting nearby agricultural areas and towns. Riverine flood has resulted in significant agricultural and livestock losses, negatively affecting local and regional economies. The study area was selected because it follows up a previous study (MMAYA & BID 2021) and the potential expansion of research regions in South America within the RECONECT project.

Figure 1(d) illustrates the study area, including nearby towns along the river. The region exhibits scattered vegetation cover and steep slopes, making it susceptible to destabilisation. The available information for this research includes a topography survey of the Arteza region and the proposed location of a protection dike according to the final report. Hydro-meteorological time series data were collected from National Meteorology and Hydrology Service of Bolivia, while data on towns, country borders, and rivers were obtained from (Geobolivia 2020).

Following the methodology of AHP, a 6 × 6 matrix was built, using the raster variables: total annual precipitation, coefficient of variation of the seasonal monthly precipitation, aridity index, slope, population density, and runoff coefficient. The matrix was computed and its statistical indexes of lambda, CI, and CR were 6.32, 0.06, and 5.8%, respectively. A CR value below 10% means the pair-wise comparison is consistent. The values of the weights obtained in the pair-wise matrix for every raster variable are 0.31 for total annual precipitation, 0.41 for variation coefficient, 0.12 for runoff coefficient, 0.08 for aridity, 0.09 for slope, and 0.21 for population density. These weights were introduced in the toolbox through a map algebra expression.

The result depends on the weights assigned to each input data obtained from the AHP. However, the weights, and hence the result, could have been different if the relative importance of a variable in the pair-wise comparison matrix were different.

Regarding the thresholds to define the four categories mentioned, they were set based on a statistical analysis from ArcGIS graphics. These statistics were extracted from the raster of the AHP, a boxplot was obtained, which the quartiles helped to define the thresholds.

The land cover computed in GEE has high accuracy, which is expressed by a Kappa training value of 99.96% and a Kappa validation value of 87.92%. The area used for croplands obtained for 2015 is 6.86 km², and for 2021 is 7.79 km² as shown in Figure 3 . This means in 2021 there was an increment in cropland areas of 11.9% compared to 2015.

The area obtained from the toolbox and ESA classification correctly identifies the brown arable land as cropland. However, the surrounding area is also cropland, but it is green with grown crops. From this it follows that ESA classification is not correct for the surrounding area because it classifies it as grassland.

The toolbox conducted spatial analysis considering soil type, land cover, erosion, maximum flow velocity, and maximum flow depth, while MMAYA & BID (2021) based their results on a hydrodynamic model. Consider as part of validation the small overlap of 30 m between the dike and the areas identified by the toolbox. It can be inferred that the toolbox successfully identified the required areas.

The main challenge with WH is how to display the final results. Raster results are tied to user satisfaction, and the criteria for this can vary among stakeholders. Some sensitivity analysis was conducted, and the results show that there are some sensitive in changing the classification of the thresholds. The obtained results are very similar results to a classification based on natural breaks. It is important to note that the underlying raster data remain constant; the only variation lies in how the information is presented to the user. Furthermore, this is a limitation to define better the thresholds to classify and display the results. Future works can incorporate diverse precipitation data, including modelled precipitation, coming from rainfall–runoff models, which use the temperature data from climate change models. This enables scenario creation and comparison against the baseline, which uses existing precipitation data from sources like measured rain gauge data, available modelled data like CHIRPS, TRMM, or PERSIANN. Another future improvement consists of using field validation data and insights from stakeholders and experts from the study area, to assign the weights of importance and compute the pairwise matrix. Finally, is necessary to refine the raster result presentation by exploring alternative classification criteria for better visualisation.

The main limitation in WR is the field validation and measurement. Therefore, the field inspection/validation is necessary not only for the results of the toolbox developed, but for any source of land cover data. This criterion is particularly relevant for Pestersko Polje, where a visit like the one made to Biesbosch National Park was not possible. Therefore, a field inspection is necessary to validate not only the results of the toolbox developed, but also any land cover data. However, we can enhance the toolbox by incorporating additional methodologies to identify various wetland degradation causes and pinpoint their locations.

The primary constraint for NRS is acquiring a high-resolution DEM from a topography survey, which serves as input for the hydrodynamic model. Nevertheless, further validation is essential through comprehensive field data collection to confirm the accuracy of flood-prone areas along the riverbank. In the coming years, these tools will be indispensable for decision-makers in managing hydro-meteorological hazards. Future iterations of the toolbox may integrate seamlessly with the hydrodynamic model, streamlining the process.

The main objective of the present research work was to enhance the existing planning toolbox by adding spatial allocation methodologies for additional large-scale NBSs to manage hydro-meteorological risk reduction. The work combines ML and spatial data processing techniques to allocate large-scale NBSs. The RF ML technique was applied to obtain high accuracy of land cover raster for WR and NRS. This approach ensures robust results for the land cover classification as defined by the statistical index Kappa, which is then used for both training and validation purposes.

To achieve this, three separate toolboxes were developed within ArcGIS to facilitate the mapping of suitable locations for NRS, WR, and WH. These toolboxes were then applied in four case study areas. Furthermore, these measures were selected based on a targeted search into the database of large-scale NBSs belonging to the EC H2020 RECONECT project, considering the preferences/needs of the stakeholders (Ruangpan et al. 2021), as well as the need to expand the previously developed toolbox.

WH is an improvement of the existing toolbox for Retention Ponds (Mubeen et al. 2021). The added methodology is based on multi-criteria analysis applying the Analytical Hierarchical Process (AHP). The pair-wise matrix computed in the process determines the weights of importance of every raster input data, meeting consistency conditions.

The WR methodology primarily relies on land cover input data. Land use change from natural areas to crops is the primary cause of wetland degradation (Yazar & Ali 2016). To identify croplands, the online remote sensing products for land cover were applied with some limitations related to spatial resolution, study area coverage (European Environment Agency 2006), accuracy, and errors in defining croplands and grasslands (European Space Agency 2020; ESRI Inc. 2021). A land cover dataset with higher accuracy was obtained using the RF classifier in the cloud platform GEE. The case studies focused on Ramsar wetlands in Biesbosch and Peštersko Polje, where areas suitable for restoration were identified. To validate the results, the toolbox outputs were compared with high-definition satellite base maps in ArcGIS. However, further verification is still needed by cross-referencing with field data.

The NRS methodology is based on multi-criteria analysis. It considers the combination of maximum velocity and maximum depth with the geomorphologic features of the study area. In this case, the HEC RAS 2D hydrodynamic model was built and simulated to obtain output rasters with peak velocities and peak depths. They were part of the inputs used in the toolbox for mapping zones for NRS. The toolbox identified destabilisation-prone areas as output results. A comparison with an existing report provided locations for a future protection dike. Decision-makers can utilise these results to determine specific areas for implementation of NRS. This approach may contribute to the protection of ecosystems, settlements, crops, and other vulnerable areas without compromising river behaviour.

The research leading to these results has been conducted within the European Union's Horizon 2020 Research and Innovation Programme under grant agreement No. 776866 for the research project RECONECT. The study reflects only the authors' views, and the European Union is not liable for any use that may be made of the information contained herein.

Regarding Bolivian case studies, the dataset used belongs to the Ministry of Environment and Water (Ministerio de Medio Ambiente y Agua MMAyA) based on the National Water Balance report (MMAYA et al. 2018), and from a technical study of preparation for pre-investment (MMAYA & BID 2021).

All relevant data are available from an online repository or repositories: (https://developers.google.com/earth-engine/datasets; https://search.asf.alaska.edu/#/; https://geo.gob.bo/; https://www.openstreetmap.org/#map=7/52.154/5.295; https://rsis.ramsar.org/; https://www.arcgis.com/apps/mapviewer/index.html?layers=cfcb7609de5f478eb7666240902d4d3d; https://land.copernicus.eu/en/products/corine-land-cover; https://worldcover2020.esa.int/viewer).

The authors declare there is no conflict.