From Forest Dynamics to Wetland Siltation in Mountainous Landscapes: A RS-Based Framework for Enhancing Erosion Control

Abstract

Human activities have caused a significant change in the function and services that ecosystems have provided to society since historical times. In mountainous landscapes, the regulation of services such as water quality or erosion control has been impacted by land use and land cover (LULC) changes, especially the loss and fragmentation of forest patches. In this work, we develop a Remote Sensing (RS)-based modelling approach to identify areas for the implementation of nature-based solutions (NBS) (i.e., natural forest conservation and restoration) that allow reducing the vulnerability of aquatic ecosystems to siltation in mountainous regions. We used time series Landsat 5TM, 7ETM+, 8OLI and Sentinel 2A/2B MSI (S2) imagery to map forest dynamics and wetland distribution in Picos de Europa National Park (Cantabrian Mountains, northern Spain). We fed RS-based models with detailed in situ information based on photo-interpretation and fieldwork completed from 2017 to 2021. We estimated a forest cover increase rate of 2 ha/year comparing current and past LULC maps against external validation data. We applied this forest gain to a scenario generator model to derive a 30-year future LULC map that defines the potential forest extent for the study area in 2049. We then modelled the distribution of wetlands to identify the areas with the greatest potential for moisture accumulation. We used an S2 mosaic and topography-derived data such as the slope and topographic wetness index (TWI), which indicate terrain water accumulation. Overall accuracy scores reached values of 86% for LULC classification and 61% for wetland mapping. At the same time, we obtained the potential erosion using the NetMap software to identify potential sediment production, transport and deposition areas. Finally, forest dynamics, wetland distribution and potential erosion were combined in a multi-criteria analysis aiming to reduce the amount of sediment reaching selected wetlands. We achieved this by identifying the most suitable locations for the conservation and restoration of natural forests on slopes and in riparian areas, which may reduce the risk of soil erosion and maximise sediment filtering, respectively. The results show a network pattern for forest management that would allow for controlling erosion effects across space and time at three levels: one, by reducing the load that originates upslope in the absence of forest cover; two, by intersecting runoff at watercourses related to sediment transport; and three, by a lack of former barriers, by trapping erosion near to the receiving wetland systems, main river axes and contributing streams. In conclusion, the proposed methodology, which could be transferred to other mountain regions, allows to optimise investment for erosion prevention and wetland conservation by using only very specific areas of the landscape for habitat management (e.g., for NBS implementation).

1. Introduction

As water is the most valuable and limited resource for biodiversity, comprehensive knowledge of water-related ecosystem functioning is essential for developing an integrated management scheme that assures its sustainability [1,2,3]. To establish priority actions on water resources worldwide, spatial assessment and environmental mapping play an essential role in defining critical areas for conservation and ecological restoration [4,5]. In the context of global change, land use and land cover (LULC) changes resulting from urban, farmland and livestock expansion [6,7] have led to the destruction of habitats, ecosystem functionality and the historical services they offer to society [8,9]. Ecosystem Services (ES) are the direct and indirect contributions that ecosystems provide to human well-being and quality of life [10,11]. They are essential for sustaining life on Earth and represent an accurate indicator of an ecosystem’s health [12,13,14,15]. The most important disturbances are related to (i) the loss of forest stands [16,17]; (ii) fragmentation of extensive, original forests in patched mosaics [18,19]; and (iii) soil erosion and water and soil contamination [20,21,22]. Therefore, understanding how ES responds to anthropogenic changes is necessary to assess and balance short-term needs with long-term sustainability goals [23,24].

Aquatic ecosystems, including wetlands and water bodies, provide different and essential functions and ES, such as habitat for wildlife, water filtration, flood protection or carbon storage [25,26]. Wetlands are typically located at defined depressions in the ground or where the water table is close to the surface [27]. There is extensive evidence that aquatic ecosystem functions can be impaired by modification of the adjacent LULC change, bringing on siltation processes at different levels [28,29]. In this regard, forests play an essential role in maintaining ecosystem functions relevant to aquatic ecosystems, such as modifying soil structure, allowing better infiltration and water retention, maintaining more stable runoff levels during drying season [30], playing a decisive role in the overall water cycle by balancing infiltration–evapotranspiration–runoff [31] and reducing soil erosion and sediment transport towards receiving streams [32].

In developed countries, woodlands and shrublands have regenerated naturally since the 19th century because of rural depopulation and land abandonment in marginal areas [33,34]. This involved the expansion of forest land and an increase in the maturity of forest ecosystems and, therefore, their capacity to provide regulating ES such as erosion protection, hydrological stability or water quality maintenance [35]. Consequently, land management strategies such as rewilding are currently at the forefront of European adaptation policies as a way to recover forest cover [36,37]. Rewilding implies passive management of ecological succession to restore natural ecosystem processes and minimise human control of landscapes. This allows for maximising environmental and socioeconomic benefits in rural areas undergoing changes such as land abandonment [38,39].

We focused on soil erosion because it is one of the biggest threats to wetland conservation, involving negative impacts such as eutrophication and increased risk of flooding and landslides [40,41]. In this regard, regulating Ecosystem Services for soil erosion prevention is paramount in mountainous areas [42]. We evaluated, in this work, soil erosion at three stages concerning natural and human drivers. The first stage consists of an initial detachment of individual particles from the soil mass, subsequent transportation by erosion agents—principally wind and water—and final siltation into land depressions or delivery to water bodies. Sediment detachment is not exclusively the result of raindrop impacts [43,44] but also originates in the overland flow as an erosive agent [45]. Both factors generate shear stress on the soil surface: if the cohesive soil strength is exceeded, sediment detachment occurs [46]. Sediment enters the river dissolved or as fragments eroded from the adjacent hills, banks or channels [47]. At the same time, streamflow provides the driving force that transports the sediment downstream within river corridors. Hence, surface and subsurface water fluxes are strongly influenced by river characteristics, including valley and channel geometry, vegetation, biota and human activities [48]. Wetland siltation occurs when the eroded material transported by river channels or overland flow settles out of the water column onto the surface as the water flow slows down [49,50]. When this occurs, the presence of sediment often significantly reduces water and the quality of the aquatic habitat for biota [51]. It makes water bodies more turbid and does not allow light to penetrate the water column [52]. Furthermore, suspended sediment can clog the gills of aquatic organisms, cover the stream bottom leading to the suffocation of fish eggs and benthic macroinvertebrates [53] and reduce dissolved oxygen levels due to higher water temperatures [54]. Additionally, pesticides, some metals and other pollutants may cling to suspended sediment in the water and increase the concentration of the contaminants and toxins [55,56]. Therefore, there is a need for direct and indirect management actions to avoid diffuse sources of sediment and pollutants, mainly related to soil erosion and runoff after rainfall events in deforested areas [57]. In cases where rainfall erosion occurs, these downstream effects are often significant and can cause severe environmental and economic damages [58].

Nowadays, nature-based solutions (NBS) are gaining considerable socio-political traction as evidence mounts regarding the multiple regulating ES and the biodiversity they provide [59]. NBS is used as an umbrella framework encompassing multiple conservation and restoration managing strategies applied over natural and semi-natural ecosystems to enhance ES provision and biodiversity [60]. The standard approach is mainly based on using natural processes and cycles of organisms, matter and energy, taking advantage of the ecosystem’s properties and following their seasonal and temporal dynamics [61]. Therefore, they are sustainable over time for potentially maintaining multiple ecosystem functions alternatively or simultaneously [62]. In this sense, NBS might be able to improve the ecosystem where they are embedded. Their effects might span across ecosystem boundaries (i.e., meta-ecosystem) [63]. As shown in [64], these properties make them particularly recommendable for managing aquatic ecosystems based on their biotic and abiotic exchanges with terrestrial habitats. To identify areas for NBS implementation, we can use different tools for semiautomatic analyses at a regional scale, characterising landscape patterns and processes at a geographical scale and temporal resolution to match the ecosystem processes of interest. In this sense, RS is a source of Earth observation data of great interest in landscape studies given its ability to follow dynamic processes [65,66]. RS has become a critical tool for landscape monitoring given its capacity to generate spatially continuous, regular and repeatable observations over large areas. It is a primary source of information for global monitoring of ecosystem processes, along with national statistics, field-based observations and numerical modelling [67,68]. Among its applications, LULC characterisation is one of the most developed branches in conjunction with geographic information systems (GIS), computational capabilities and statistical techniques developed over the last decades. All have turned ecological models into powerful tools for spatial analysis aimed at (i) exploring the various mechanisms that force changes in LULC, as well as their social and economic driving forces, (ii) projecting the potential environmental and socioeconomic impacts derived from these changes across future scenarios, and (iii) assessing the influence of policy alternatives and management regimes on ecosystems functionality [65].

Much of the work to protect aquatic ecosystems, particularly from siltation effects, is based on creating vegetated buffer strips for stream management [69,70]. Beyond this, the quality of the services generated depends largely on the conservation state of forest stands [71,72]. Therefore, the ecosystem function is defined by the extension, structure and composition of forest patches and the optimal distribution across the landscape of these vegetation barriers [73,74,75]. However, very few publications address this problem from a comprehensive point of view, including reducing sediment production in the headwaters, sediment transport in river channels or diffuse runoff and final siltation in receiving water bodies.

2. Objectives

The main objective of this work is to apply an RS-based framework to identify suitable locations related to the conservation and restoration of natural forests in hillslopes and riparian areas that reduce the risk of soil erosion and increase sediment filtering, respectively. This approach will allow analysis of the vulnerability of siltation on wetland ecosystems in current and future LULC scenarios concerning soil erosion, sediment transport and sedimentation. To do so, we carried out a series of specific steps:

• Mapping LULC dynamics from 1989 to 2019 by coupling Landsat and Sentinel 2 imagery.
• Modelling forest distribution in a future scenario for 2049 based on a reliable maximum forest cover increment defined for the study area.
• Mapping current wetlands distribution supported by extensive fieldwork and photo-interpretation.
• Modelling landscape processes related to the production, transport and deposition of sediment.
• Application of multi-criteria analyses to identify potential areas for the implementation of NBS (i.e., forest restoration and conservation) to reduce erosion effects on wetland ecosystems.

3. Study Area

The study area covers the Sella, Cares and Deva catchments, which extend across the northern slopes of the Cantabrian Cordillera (Spain) with a total surface area of 2485.20 km² that spans more than 400 km parallel to the Cantabrian Sea [33] (Figure 1). The current morphology results from small-scale tectonic movements and the action of wind, fluvial and glacier erosive agents, thus forming very active geosystems that gave rise to a great diversity of landscapes [76]. The climate is humid, oceanic and temperate [73], and the region is subdivided into different biogeographical provinces distributed along the Colino, Montano, Subalpine and Alpine bioclimatic floors [77]. Combined with steep terrain and sharp altitudinal gradients, they make the mountain range a longitudinal, and sometimes a transversal, border for species transit [76]. This heterogeneity of conditions means that the region is home to one of the greatest botanical richnesses of the European mountain regions [78], as well as an enormous wealth of fauna, including 18 species of amphibians, 22 reptiles, 190 breeding birds and 67 mammals [79,80,81]. In addition, the above-mentioned border effect translates into a high percentage of endemism [82]. Due to the borderline nature, i.e., where different ecosystems converge between the Atlantic and the Mediterranean regions, the area is highly sensitive to climate change [76].

Short and steep rivers with great erosive power [83] and frequent floods of varying magnitude [84] characterise the area’s hydrology. The biodiversity of aquatic systems in the region, including river networks and lentic elements such as lagoons, lakes, springs and peat bogs (i.e., wetlands from now on), show a high richness of diatoms, macroinvertebrates and amphibians [76,78]. Local endorheic lakes formed by the action of glacial ice during the Pleistocene [85] created complex hydrological networks with streams, wetlands and other water bodies which are highly connected and interdependent [86]. Wetlands are more complex to identify and characterise, and are usually former glacial lakes or the result of an outcrop of water in flat terrain, which is very common due to the predominance of limestone materials throughout the territory [26,86].

The land cover gradient is characterised by the legacy of human management and land-use practices, which dates from the Neolithic [87] and whose influence on the forest has been prominent since the 13th and 14th centuries [88]. In the last 400 years, this influence has been based on deforestation for stockbreeding, agriculture and wood production through the combined use of fire and traditional agricultural techniques. These heavily managed landscapes [88,89] suffered further transformations during recent decades because of land abandonment and rural exodus [90], with the subsequent increase in the occupied area and density of forests and shrublands [91,92] and the decrease of crops, pastures and grasslands [89,93]. In contrast, large parts of the study area have not experienced major historical deforestation and present mature Eurosiberian forest patches of beech (Fagus sylvatica), birches (Betula ssp.), oaks (Quercus robur) and Mediterranean vegetation in lowlands (Quercus pyrenaica or Q. ilex). Here, the presence of Brown Bear (Ursus arctos) and Cantabrian Capercaillie (Tetrao urogallus cantabricus) demonstrates a better landscape conservation status [76,94,95].

The combination of these factors related to relief, hydrology and LULC change has made the numerous aquatic ecosystems in the study area vulnerable to siltation. The on-site impacts of soil erosion are widespread and manifest in very high rates of soil loss, in some cases exceeding 200 T·ha⁻¹·y⁻¹ (INES), generating off-site costs like siltation of dams and ponds and downslope damage to property, roads and other infrastructure [96,97].

4. Materials and Methods

We developed a flow diagram to summarise the main methodological steps connecting forest dynamics to wetland siltation in mountainous landscapes (Figure 2).

4.1. In Situ Data across Environmental Gradients

To train LULC and wetland ecosystems classifications, we obtained a data set of reference points distributed across the study area. Those included field campaigns carried out from 2017 to 2021 and photo-interpretation using high-resolution aerial photographs from the PNOA program (20 cm of pixel size), freely available for the whole of Spain for the same years (CNIG; https://centrodedescargas.cnig.es). We paid particular attention to selecting accurate and representative reference data across the whole territory, avoiding sampling noise from border-type effects and mixed categories, as well as reducing sampling bias about location, slope and orientation [98]. Some data were also validated through Google Street View to assure this [99,100]. We achieved this by collecting training points for LULC classification by using a 2 km grid to stratify the area, reducing spatial autocorrelation and pseudoreplication in medium-resolution images of heterogeneous mountain areas [101]. We accounted for 1.821 field-checked locations: 457 for broadleaf and 205 for coniferous forest, 320 shrublands, 227 pastures and around 700 for other land cover types. On the other hand, we collected 282 reference points (GPS quality) of wetlands across the study area, with varying soil moisture gradients and topographic constraints. Checking these reference locations, we also digitalised wetland patches visible during low-level flows by using high-resolution aerial photographs from the end of summer, getting the reliable size and shape of 141 wetlands for subsequent distribution modelling. Absence points were selected randomly across the area following three criteria: (i) excluding the former 141 photo-interpreted patches, (ii) only within shrublands and grasslands obtained in the previous LULC classification and (iii) with a minimum distance of 200 metres among them to avoid pseudoreplication.

4.2. Satellite Data and Pre-Processing

LULC mapping was performed by applying a supervised classification over a cloud-free satellite-imagery mosaic covering the whole region, using different scenes from July 2019 to assure a maximum sun elevation angle and less cloud coverage. Top of the atmosphere (Level-1) data from the Copernicus Sentinel-2 constellation were chosen for this task. Data were acquired by two operational satellites (S2A and S2B) that provide multispectral images at 10, 20 and 60 m of spatial resolution, at a five-day revisiting time. Copernicus S2A was launched in June 2015, followed by its twin satellite (S2B) in March 2017. The images were downloaded from the Copernicus Open Access Hub service (https://scihub.copernicus.eu/, accessed on 21 September 2020) and atmospherically and radiometrically corrected to the Bottom of Atmosphere (BOA) level using the Sen2Cor algorithm available in the SNAP Sentinel toolbox. Using the bilinear interpolation algorithm, all image bands were subsampled to a 10 m grid. The subsampled bands (b5, b6, b7, b8A, b11 and b12) were then stacked with original 10 m bands (b2, b3, b4 and b8) for deriving spectral indices related to vegetation functional traits such as productivity and phenology. In particular, we obtained the Normalised Difference Vegetation Index (NDVI), commonly used in land cover mapping for enhancing vegetation productivity and biomass [102,103]. Some studies report that LULC classification benefits from incorporating information from the red edge and the near-infrared of the electromagnetic spectrum by exploiting S2 bands 4 and 8. The Normalized Difference Green Vegetation (GNDVI) index also exploits the visible green spectrum (540–570 nm). Apart from being more sensitive to chlorophyll concentration than NDVI, it is good senescence and leaf hydric stress indicator [104] that occurs during the summer period in the study area. A complimentary Digital Elevation Model (DEM) of 5 m of spatial resolution was obtained from LiDAR data, also available for the whole country at CNIG (https://centrodedescargas.cnig.es), and resampled to 10 m and 30 m to match S2 and Landsat pixel sizes, respectively. In turn, LULC classification in the past was developed over a time series of 41 cloud-free Landsat images collected from 1989 to 2019, processed with NASA’s Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS). Such data were obtained from the US Geological Survey (USGS) Earth Explorer archive (http://earthexplorer.usgs.gov/, accessed on 21 September 2020). Similar to the previous process, a minimum amount of cloud cover and high elevation angles were required, which focused the search on images taken during the summer period of each year.

We used the same S2 mosaic to map wetland distribution, but corrected it with the iCOR algorithm as it presents a better performance for water-related ecosystems [105,106]. We generated two specific spectral indices related to water presence—the Normal Difference Water Index (NDMI) and the Automated Water Extraction Index (AWEInsh) [107,108,109]. Furthermore, we created the Topographic Wetness Index (TWI), derived from topographical attributes, to define the areas with the highest potential for moisture accumulation [110].

4.3. Mapping LULC Trends and Wetlands Distribution

LULC classification has been made for three time periods spanning 30 years (1989–2019–2049). The beginning of any supervised classification involves the definition of the legend, i.e., in which categories the image is going to be classified.

Table 1. Definition of basic legend for LULC classification [91].

Legend		Description
1	Broadleaf forest (BLF)	Vegetation formation composed mainly of mature trees with shrub and undergrowth presence in which broadleaved species predominate. Applicable to mature forests of natural or anthropogenic origin such as pure or mixed broadleaved stands; riparian and gallery forests; broadleaved deciduous thermophiles forests; evergreen sclerophyllous forests; arborescent thickets with sclerophyllous broadleaved species; logging areas; young broadleaved tree plantations; also, to young broadleaved plantations.
2	Coniferous forest (CNF)	Vegetation formation consisting of mature trees with shrub and undergrowth presence in which coniferous species predominate. Applicable to mature forests of natural or anthropogenic origin such as pure or mixed coniferous stands; deciduous, coniferous forest; arborescent scrub. Furthermore, to young conifer plantations and coniferous wooded dunes.
3	Shrublands (SSH)	Areas of dense, herbaceous vegetation with occasional scattered trees. Representative of natural development of broadleaved and/or coniferous forest formations in transition by natural succession on abandoned agricultural land; forest regeneration after damage (fire, landslides); stages of forest degeneration caused by natural or anthropogenic stress factors; reforestation after logging or planting in previously non-forested natural or semi-natural areas.
4	Agricultural land, crops (AGR)	Areas of cultivated land in agricultural use for arable crops irrigated periodically or permanently, using permanent infrastructure (irrigation channels, drainage network and additional irrigation facilities). Most of these crops cannot be cultivated without an artificial water supply.
5	Pasture and grassland (GRA)	Permanent grasslands characterised by agricultural use or heavy human disturbance. They have a floral composition influenced by human activity, usually grazing or mechanical harvesting of grass meadows.
6	Denude rock and bare soil (DEN)	Areas of boulders, cliffs, rocky outcrops, including areas of active erosion, rocks and natural no vegetated expanses of sand or pebble/gravel. In coastal or inland locations, such as beaches, dunes, gravel, including torrential stream channel beds. It also includes areas with sparse vegetation consisting of herbaceous and/or semi-woody species covering 10–50% of the area. Includes steppes, tundra, lichen heaths, wastelands, karst areas and sparse vegetation at high altitudes.
7	Urban and human-derived areas (UHD)	Areas of continuous or discontinuous urban land, as well as industrial areas, road networks, port areas and artificial areas occupied by extractive activities. Areas with artificial vegetation are also included.
8	Water ecosystem (WAE)	Inland sea area and lower tidal limit as well as lakes, ponds and flowing waters formed by all rivers and streams, including also artificial water bodies such as reservoirs and canals.

shows the legend for the present study based on the typologies used in [111] according to vegetation physiognomic units that assure an optimal classification with medium resolution sensors in heterogeneous and mountainous environments.

We then applied a per-pixel supervised classification using a Support Vector Machine algorithm (SVM) which collected in situ data and the predictors explained above. The SVM algorithm generates a hard map, assuming homogeneous pixels across the landscape continuum and creating a complete classification for the S2 mosaic. When one pixel does not belong entirely to any defined class, the algorithm assigns it as “not evaluated” or to the category with the highest probability of belonging [112]. Therefore, the process calculates an uncertainty that can be evaluated by complementary soft-classification approaches and used to refine LULC patterns and trends [111].

The analysis of LULC dynamics requires assuring spectral coherence between Landsat and Sentinel images. Radiometric differences must be considered due to aspects such as phenology and vegetation growth rates, varying atmospheric conditions, soil moisture or the presence of aerosols that may modify BOA reflectance and therefore affect the post-classification comparison of LULC patterns. Without a high temporal resolution of LULC maps over time, the performance of an adequate atmospheric and topographic correction, the quality of the training data or even the selection of images during the summer period do not guarantee LULC comparability over time, which may result in inconsistent change detection [65,101]. Here, we applied a sensitivity analysis that calculates the averaged values of NDVI for each LULC category defined in the legend for the year 2019 and searched the most spectrally similar image in the past, even when dealing with different sensors, thus guaranteeing comparability of maps over time. We then calculated the NDVI value for each in situ data point for the whole series of Landsat and S2 images, together with the 25th and 75th percentiles between which the interquartile range is found, i.e., 50% of the sample values (see Figure 3). These averages are then calculated for each category and Landsat image in the time series, identifying all in situ data points within the current image’s interquartile range. In this way, the June 1989 scene was selected as the one that assures a larger consistency in the time series for the LULC change analysis. After selecting the scene, LULC classification was performed similarly to S2, applying the SVM algorithm and using only the reference points of each category with NDVI values within the interquartile range to keep the spectral unmixing of LULC categories as much as possible. Current S2 maps and ancillary topography data were finally rescaled to 30 m using a nearest neighbour approach to match the Landsat data.

After classifying past and current LULC patterns, we generated a matrix of change by comparing the maps from 1989 and 2019. We established an increment rate of the forest cover and applied the Scenario Generator Proximity Based model of InVEST (Integrated Valuation of Ecosystem Services and Trade-offs, see https://naturalcapitalproject.stanford.edu, accessed on 17 May 2021) [113,114,115] to create a future BAU (Business As Usual) scenario allocating the maximum increase in the forest category. Previously, we merged BLF and CNF typologies (see Table 1) into a single forest category called FOR to avoid the uncertainty generated by the translocation of pixels of both typologies without any specific analysis of the phenology, drivers and abiotic limiting factors at the species level. The InVEST service shows that this tool requires establishing the focal LC that sets the proximity rules from which the scenario will be determined—in this case, SSH and GRA. The scenario generator will then convert pixels of grasslands and shrublands into forests, starting from the edge of focal patches at the base of the higher probability to be modified [113,114].

We used the same supervised classification with an SVM algorithm for mapping wetland distribution over a combination of the S2 imagery corrected by the iCOR package, spectral indices and ancillary topographic variables [116]. The model output provides per-pixel probabilities of wetland occurrence with continuous values ranging from 0 to 1. With the previously digitalised wetland patches, we extracted a suitable threshold for creating presence–absence maps from model outputs and used it as a cut-off threshold that allowed us spatially define wetland areas.

4.4. Potential Erosion, Surface Runoff and Deposition Areas

The potential erosion in hillsides and sediment delivery on wetlands ecosystems was determined using the NetMap tool with an automated process that analyses watershed functioning according to characteristic parameters such as topography and pluviometry [117,118]. For its use, we generated a Synthetic River Network (SRN) from the digital elevation model (DEM) with 10 m of spatial resolution to match the S2 data. This network is an optimal hierarchical spatial framework to articulate hydrological and environmental information [119]. The stream networks are delineated using flow directions inferred from the DEM and the algorithms described by [120]. We applied drainage enforcement in lower relief areas (slope less than 30%) by reducing the elevation of the current cells in the DEM by two metres, using the current locations of reference river channels to avoid these cells acting as sinks. Nodes divide the network where several river channels meet, as these junction areas imply changes in channel morphologies [117,118]. Then, we applied specific NetMap packages that encompass shallow landslide, gully and surface erosion processes through geomorphologic criteria to obtain the Generic Erosion Potential index (GEP). GEP provides a relative measure of potential erosion at the source (e.g., from 0 to 1) calculated from the flow accumulation to each pixel [118,121]. The GEP index will give larger values in steeper and convergent slopes [118], as they are more prone to landslides and surface erosion [122]. In addition, the RUSLE Length–Slope (LS) factor was added as a measure of the sediment transport capacity of surface runoff [123]. Finally, we defined deposition areas as plains with a slope of less than 5% [124,125,126] that intersect with the river network which coincides with the wetlands defined in the previous classification.

4.5. Multi-Criteria Analysis: From Forest Dynamics to Wetland Siltation

We designed a multi-criteria analysis to identify priority areas to reduce soil erosion on hillsides and the probability that the sediment generated will be transported and finally deposited in the wetlands of the study area. This action targets forest conservation and restoration actions, but specifically also spontaneous vegetation regeneration following secondary succession processes across hillslopes and riparian areas (i.e., valley bottoms and wetlands receiving runoff from adjacent zones). To this end, we analysed each process separately.

For evaluating erosion at the source through the GEP, we selected the surfaces with a value greater than 0.50 upstream of the wetland ecosystem of interest [121,127] with an influence area of 30 m to the transport channels [127] (Figure 4). We selected all upstream channels within the drainage area for transport issues, discarding the rest. All overland flow transport channels obtained through the LS index have been included; we identified those pixels with an LS value greater than 0 [123]. We then define four categories of transport capacity based on the quantiles of the LS index, with values from 1 to 74 (low, medium, high and very high), where the channels of the river network are always categorised as “very high”. The influence area of the sediment transport network is 5 m for all those with a slope greater than 30% and 10 m for the rest (Figure 4) [127,128]. Finally, the sedimentation process is allocated to the entire area within 20 m of the target wetlands with a slope of less than 5% (Figure 4).

According to [64], the spatially explicit delineation of erosion at source, sediment transport and sedimentation processes constitute the area of abiotic potential flow. The presence of forest systems in this area regulates the transfer of sediment between the terrestrial environment and the wetlands, so different areas for forest management can be established, either at the maximum (i.e., covering all possible erosion processes) or only in those areas where the problems are more important (see results section). In this sense, the area of abiotic potential flow was overlaid with: (i) the current vegetation map (2019) to identify forest areas that are currently preventing erosion, and (ii) the future forest expansion map (2049) to identify the areas that may experience a reduction in erosion due to future LULC dynamics.

5. Results

We obtained the distribution and dynamics of LULC typologies according to the legend of Table 1, which show a significant increase in forests over time, from a coverage of 14.84% in 1989 to more than 38% in 2019 (see

Table 2. LULC distributional areas for the past, current and future scenarios.

L.U.		1989				2019				2049
		Surface (ha)		% Surface		Surface (ha)		% Surface		Surface (ha)	% Surface
BLF	FOR	36,499.40	36,874.69	14.69	14.84	94,517.00	96,499.57	38.04	38.84	153,799.57	61.89
CNF		375.29		0.15		1982.57		0.80
SSH		99,809.43		40.16		85,214.25		34.29		39,627.61	15.95
AGR		21,619.85		8.70		11,090.19		4.46		11,090.19	4.46
GRA		57,229.07		23.03		21,533.37		8.67		9820.01	3.95
DEN		31,819.08		12.80		31,449.75		12.66		31,449.75	12.66
UHD		767.20		0.31		2332.19		0.94		2332.19	0.94
WAE		380.01		0.15		380.01		0.15		380.01	0.15

and Figure 5). This represents a total area of 59,624.88 ha, i.e., 1987.50 ha/year, which was used to model forest expansion using the InVEST tool. The overall accuracy for the current scene is 86.4% ± 1.4 with a 95% confidence interval. The classification was able to identify similar LULC typologies, getting, e.g., minimum Jeffries–Matusita separability values of 1.15 between the AGR and GRA coverages, and 1.56 SSH and GRA. In the 1980s’ classification, we obtained an 80.2% ± 2.3 overall accuracy with a 95% confidence interval.

From 2019 to 2049, the forested area is predicted to increase by more than 57,000 ha, which means shifting from occupying less than 40% to more than 61% of the total area. This large increase in the forest category for 2049 is related to the decrease in shrublands. The process has two stages: first, a slight reduction between 1989 and 2019, from 40% to 35% of their total area or about 15,000 ha in total or 480 ha/year, and secondly, the extreme reduction simulated between 2019 and 2049, from 35% of the area to less than 16%, with about 45,000 ha or more than 1500 ha/year. This effect is offset by pastures and grassland, with a significant reduction from 1989 to 2019 of more than 36,000 ha, or 1200 ha/year, and a slight decrease between 2019 and 2049 of around 20,000 ha, or 660 ha/year. In addition, agricultural land has a reduction between 1989 and 2019 of about 10,000 ha and then remains constant in the model until 2049. Urban areas have increased from 1989 to 2019 and remain stable until 2049. Finally, the water ecosystem and denude rock/bare soil remain constant throughout time.

The wetland modelling outcomes were able to identify 171 reference points across the study area from the 282 available, leading to a detection capability of 60.6% ± 6.8 with a 95% confidence interval. In turn, the 141 wetlands mapped by expert photo-interpretation occupied a total area of 170.50 ha. Wetland models, after applying suitable thresholds over continuous 0–1 probability values and acknowledging the capabilities of such an RS-based modelling approach, occupy a total area of 395.66 ha. This results in an overlay percentage of 38.54% or, in other words, an overestimation of wetland models in the expert criteria.

The results obtained for GEP indicate a greater erosion potential on steeper slopes, which tend to coincide with non-forested headwater basins. In these areas, sediment flow channels show that temporary water networks or runoff flows transport the most considerable amount of sediment towards the wetland ecosystems as opposed to the permanent networks, which make a much smaller yet significant contribution. The NetMap tool also allowed identifying which reaches receive the most important contribution of sediment (see Figure 4).

Finally, the multi-criteria analysis identified optimal areas to apply forest management strategies based on conservation and passive restoration actions (692.13 and 1289.11 ha, respectively). These current and future forest patches may allow for reducing erosion process at three stages (

Table 3. Area available for the implementation of forest management actions (i.e., NBS of forest conservation and forest restoration) in relation to the analysis of erosion process.

	Erosion at Source	Sediment Transport				Sediment Deposition	Total
		Low	Medium	High	Very High
Forest conservation (ha)	92.48	308.22	122.51	85.77	77.17	5.98	692.13
Forest restoration (ha)	137.48	556.24	185.90	108.47	236.42	64.6	1289.14

): (i) reducing erosion production at the source through the presence of forests on drainage slopes, (ii) filtering sediment inputs and retaining sediment transport in permanent and temporary water networks using riparian forest strips, and (iii) filtering out sediment delivery or deposition to wetlands through riparian forest buffers in surrounding areas. The larger areas to preserve and restore correspond to sediment transport (see Table 3). Dealing with the intensity of sediment transport, larger areas correspond to high transport capacity, or in other words, to those where surface runoff is more important. Secondly, erosion at the source, to which forest management actions could be applied, affects less than 15% of the area. Finally, the areas addressing sediment deposition occupy 3.56% of the area of potential erosion flow. As a graphical example, Figure 6 shows potential areas for establishing forest management actions in one of the studied wetlands.

6. Discussion

In this work, we applied a novel RS-based modelling approach that couples forest dynamics and potential wetlands siltation to identify optimal areas for implementing NBS that regulate sediment flows in mountainous systems.

6.1. Landscape Modelling Approach

LULC mapping has been based on sensitivity analysis to detect spectrally comparable satellite images from different sensors within a time series of Landsat and S2 scenes. Various studies showed that both sensors are directly comparable after a TOA correction when the images were taken on the same date [129,130]. Other studies developed a Time-Series-Based Reflectance Adjustment (TRA) approach for reducing the reflectance differences between both sensor observations after a TOA correction [131] and linear regression models to adjust reflectance differences between sensors for each pixel and spectral band. We developed here a simplified, conservative approach that maintains the original values of band images. To achieve this, we based our analyses on a sensitivity selection by NDVI according to the LULC legend to obtain comparable data at the BOA level.

The RS-based modelling of wetlands allows assessing the uncertainty of boundaries extracted by automated algorithms in complex landscapes against expert-based data collected on the ground [101,132,133]. These RS-based models’ uncertainties can be linked to several issues. One deals with the nature of wetland typologies. Among the 282 wetlands reference points, we have different typologies including lagoons, lakes and peat bogs. With the 10 m pixel size of the S2 images, we are, in principle, able to identify some of these aquatic systems [134,135]. Nevertheless, other aquatic ecosystems such as cattle watering holes or springs are locations where the image pixel’s size and spectral uncertainty may preclude mapping their occurrence on the ground. A second reason can be related to the need for the temporal dynamics (i.e., phenology or seasonal trends) of satellite information to catch up with the systems’ properties that allow for reliable mapping exercises. In [136,137], among studies, has been demonstrated the capability of using satellite image time series to map and monitor biodiversity patterns at a large scale. Moving beyond the validation of wetland occurrences (i.e., made here by overlying modelled patterns and reference points), we assessed the method’s capability for matching the spatial pattern of wetlands mapped using expert criteria and aerial photographs. In general, we observed an overestimation of the models over photointerpretation patterns. Nevertheless, assessing local boundary uncertainty or fuzzy boundary width in wetlands may result in artificial and not wholly reliable in the base of existing images of known class characteristics and locations, and this may hamper discrimination and variability [138,139]. Discrimination consists of the ability of the photointerpreter (or an algorithm) to detect a difference in texture. Variability consists on the intrinsic spatial variability of the texture itself because of the fuzzy boundaries [138,140]. Acknowledging these complexities, we determine that RS-based modelling of wetlands may help examine the effect of context on boundary uncertainty in complex landscapes, i.e., how one might assess the uncertainty of boundaries extracted by automated algorithms based on RS against reliable data collected on the ground to advance more reliable and automatic mapping procedures. Nevertheless, this issue needs further research to determine the omission and commission errors of both techniques.

The combination of LULC and wetland distributions in this work allowed assessing the potential flow of sediment following a meta-ecosystem perspective [63]. This identified, on the one hand, all potential pathways of sediment through terrestrial landforms to wetlands, but also depicted where terrestrial ecosystems are connected to these abiotic flows in order to reduce the magnitude of the sedimentation processes in target wetlands. In this sense, erosion at the source, sediment transport and sediment deposition has been analysed independently. This is a major advance over other approaches that either consider only some of the processes in isolation, or cannot spatially discretise them to assess accurately in which terrestrial structures they occur [141,142,143]. In this way, our methodology identifies not only different forest habitats (e.g., riparian forests or hillside forests), but also different managing strategies (e.g., conservation or restoration), that could be applied at various scales and locations to reduce the sediment transfer between terrestrial and wetland ecosystems. Therefore, using the proposed approach, we could first prioritise which locations need urgent actions and which are not at risk from erosive sediment inputs. Furthermore, we could identify over which process and locations we can act more effectively, e.g., erosion production at the source, sediment transport by rills and channels or sedimentation in flat areas.

6.2. Implications for Landscape Management

Our results show that, in the study area, most of the sediment is generated in the headwaters. This is in agreement with other studies carried out in mountainous territories where slope steepness and convergence cause surface and subsurface flow to become concentrated and thus contribute to soil erosion [117,118,144]. In this sense, vegetation is an important agent for stabilizing steep regolith-mantled slopes. Vegetation reinforces the regolith with roots [145,146] and, to a lesser extent, modifies soil moisture and subsurface hydrology through transpiration, canopy interception, the redistribution of rainwater and the development of preferential flow paths via root systems [147,148]. From an ecological point of view, more complex ecosystems and mature vegetation contribute to greater protection of the soil properties [40]. In fact, woody vegetation stabilizes soil surfaces by impeding soil movement, generating lower losses in comparison to other vegetation types [149,150,151]. In our study area, forest conservation actions would not affect large extents (around 4% of the total area implied in sediment flows), so fully preserving its current adequate condition should be a priority. However, as has been quantified by [64] for the same study area, an important part of the sediment is generated in areas where forests are not present, now nor in the future, because of limiting abiotic factors such as altitude, slopes and climatic conditions [101]. In this sense, the total area of sediment generation is significantly larger than the 229.96 ha (Table 3) identified in our study as a potential area to reduce the production of sediment in wetlands. In these cases, management should focus on reducing sediment transport and deposition by means of, e.g., riparian vegetation, in order to regulate sediment flows that are difficult to mitigate at the source. For example, according to some studies, riparian areas covered by dense and mature forests were ≥75% effective in trapping sediment before reaching a water body in agricultural land [127], and 53% to 96% in piedmont areas [152,153]. For this purpose, our methodology allows us to locate not only the temporary and secondary channels with a high capacity to transport sediment from the production areas to the wetlands, but also the areas where sediment arrives at the wetlands (Table 3). However, during high-intensity precipitation events, which are expected to increase in frequency because of climate change [154], the effectiveness of these buffering systems may be limited by the concentration of the surface flow of water and sediment before they even reach the delivery channels [155,156,157]. Under these circumstances, conservation and restoration of hillside forest in areas of low and medium transport capacity become of paramount importance (Table 3): the increase in the roughness of the initialization rills of water and sediment flow can contribute to making it more diffuse, which increases the effectiveness of the riparian filters [158].

In terms of landscape dynamics, we observed in the study area a widespread and intense forest expansion that can be read as an opportunity to provide erosion regulation [159,160]. Our results about forest dynamics were compared to another study in a similar area [87], which provides an average annual forest cover increment rate for the period between 1956 and 2004 of 0.54%, representing an increment of 1.94 ha/year. If we extrapolate this value to our study area, we get a forest expansion rate of 1342.01 ha/year. With these data, we estimate that the average area occupied by forests for the year 1983 was 32.01% and 47.53% in 2004, which implies an increment of 2.94 ha/year. Extrapolating these rates to our study area, it represents 2029.58 ha/year. This reinforced our estimation of future LULC for the year 2049, which is also coherent with the dynamics observed in other mountain regions [111,161,162]. In the base on these results, the reactivation of secondary succession processes would allow creating, in the next 30 years, a total of 597.01 ha of new forest stands for the study area, which represents a 115% increase compared to the current forest that, in turn, is linked to the erosion regulation and subsequent sediment delivery to the wetlands (see Table 3). Considering these trends, therefore, passive management of forest expansion (i.e., rewilding) seems to be a very cost-effective adaptive strategy for improving wetland conditions and functioning. Apart of that, other restoration techniques such as sediment trapping [163] or human-made afforestation could be also applied when erosion effects are more important. This would assure wetland conditions in case a wetland is close to clogging or receives sediment with a large amount of nutrients [164], or simply when future forests will not be able to reach specific critical locations at the base of our potential simulation [91]. However, several authors pointed out the negative effects that rewilding may have on preserving the rich diversity of European cultural landscapes [165,166], as well as its impact on traditional agrarian activities and, consequently, some food-provisioning ES [167]. In fact, in our study area, forest expansion is driven by a decrease in livestock pressure, which may involve the loss of highly productive grasslands [94]. In this context, our methodology could contribute to decision-making by identifying areas where the process of secondary succession is most favorable to wetland conservation with minimum area requirements to be managed by NBS implementation (e.g., rewilding). Furthermore, the framework is suitable to be inserted in integrative approaches that seek to design effective landscapes in response to not only social demands but also environmental conditions of the territory (e.g., Blue and Green Infrastructure Networks) [64]. In this sense, the erosion regulation ES for wetland protection would be considered together with other ES such as flood regulation, food provision or aesthetic values. By simultaneously assessing factors associated with wetland conservation such as biodiversity or water body connectivity, and the relationships with other ES (e.g., synergies or trade-offs; [167]), the areas of the territory for NBS implementation would be prioritized in a multivariate and cost-effective manner.

7. Conclusions

The combination of information about LULC dynamics, wetland distribution and erosion processes has allowed establishing an innovative spatially explicit RS-based workflow that allows addressing potential ecological and hydrological problems of wetlands in mountainous environments by using NBS related to forest ecosystems. In the context of this work, ecosystem management actions should be focused on conserving current forest stands and applying passive restoration techniques in areas where future forests could be established near receiving areas where the delivery of sediment from streamflow or surface runoff may reach the aquatic ecosystems of interest. The general purpose of the stratified approach would tackle not only slowing down erosion processes in the short to long term, but doing so more efficiently, thus ensuring a better ecosystem functionality that could be executed under limited available funding opportunities for ecological restoration.