Modelling Sediment Retention Services and Soil Erosion Changes in Portugal: A Spatio-Temporal Approach

Abstract

Soils provide important regulating ecosystem services and have crucial implications for human well-being and environmental conservation. However, soil degradation and particularly soil erosion jeopardize the maintenance and existence of these services. This study explores the spatio–temporal relationships of soil erosion to understand the distribution patterns of sediment retention services in mainland Portugal. Based on Corine Land Cover maps from 1990 to 2018, the InVEST Sediment Delivery Ratio (SDR) model was used to evaluate the influence of sediment dynamics for soil and water conservation. Spatial differences in the sediment retention levels were observed within the NUTS III boundaries, showing which areas are more vulnerable to soil erosion processes. Results indicated that the Region of Leiria, Douro and the coastal regions have decreased importantly in sediment retention capacity over the years. However, in most of the territory (77.52%), changes in sediment retention were little or were not important (i.e., less than 5%). The statistical validation of the model proved the consistency of the results, demonstrating that the InVEST SDR model is an appropriate tool for estimating soil loss potential by water at regional/national levels, although having its limitations. These findings can be relevant to support strategies for more efficient land-use planning regarding soil erosion mitigation practices and to stimulate further investigation at a national level on this important ecosystem service.

1. Introduction

Soil erosion is a natural process responsible for shaping the physical landscape through the distribution of weathered materials produced by geomorphic processes [1]. However, when soil erosion occurs in an accelerated rate due to anthropogenic activities, wind or water, deterioration or loss of the natural soil functions is likely to ensue [1]. Soils perform a range of key functions, including food production, storage of organic matter, water and nutrients cycling, and habitat quality for a huge variety of organisms [2]. Preserving soil resources through erosion prevention is a safeguard procedure to protect the ecological environment and the ability of soils to contribute to ecosystem functioning [3].

Soil loss by water is closely related to rainfall, snowmelt and irrigation, partly through the detaching power of water drops striking the soil surface, and, partly, through the contribution of rain to runoff [2]. Soil erosion by water has become one of the greatest global threats to the environment [4]. As a consequence, soil condition, water quality, species habitats and the provision of ecosystem services are negatively affected, which highlights the importance to quantify the impacts of soil erosion by water and developing effective measures for soil and water conservation [5]. Due to the difficulty to measure soil erosion at large scales, soil erosion models are suitable tools for regional and national estimates [6]. However, the high heterogeneity of soil erosion causal factors combined with often poor data availability remains an obstacle for applied conservation strategies [6].

Using a combination of remote sensing, Geographic Information Systems (GIS) modelling and census data, several studies have demonstrated the effects of land use and land cover on soil erosion globally [3,7,8,9,10]. In Europe, one study explored the use of the European Soil Erosion Model (EUROSEM) to simulate erosion processes, explicitly for rill and inter-rill flow [11]. Also, the RUSLE2015 model estimated soil loss at 100 m resolution for the European continent [12]. In Africa, a recent study analysed soil loss and sediment exportation at the Winike watershed in Ethiopia, concluding that land-use changes greatly affect the amount of soil loss in cultivated areas [13]. In China, authors evaluated the soil erosion at a regional scale at Yunnan Province using the Chinese Soil Loss Equation (CSLE), which allowed a more accurate soil erosion map for that province [14].

Particularly for Portugal, some studies have been carried out for modelling soil erosion at local scales (e.g., [15,16,17]). One study was about nutrient retention by trade-offs between sediments and vegetation types in Ria de Aveiro lagoon (central Portugal) [15]. Another investigation analysed the effects of land abandonment on soil erosion and land degradation in the River Côa Valley (north-eastern Portugal) [16]. Recently, a study determined the influences of gully erosion in steep regions in the northern territory of Portugal [17]. Albeit these studies have been made in different regions of Portugal, a deeper and validated study is yet to be carried out to explain the effect of sediment retention on soil erosion in the entire territory. To fill this gap, the overarching goal of this research is to explore the spatio–temporal distribution of soil erosion by understanding the spatial patterns of the sediment retention capacity in mainland Portugal, based on Land Cover changes from 1990 to 2018. Specifically, it aims to: (i) estimate the soil loss at a pixel scale, and to (ii) estimate sediment retention variations at NUTS III level.

Using the InVEST Sediment Delivery Ratio (SDR) model to determine the behaviour of sediment retention in Portugal’s mainland, the results provide a unique perspective on soil erosion and sediment retention for Portugal, contributing useful information to design effective landscape planning for soil and water conservation.

2. Materials and Methods

2.1. Study Area

This study focuses on mainland Portugal (Figure 1). Portugal is a country in southern Europe, occupying a total area of 92,212 km², whereas the mainland has a total area of 89,102.14 km², with 23 statistical boundaries defined as NUTS III [18,19]. The mainland is located on the southwest of the Iberian Peninsula, bordering Spain to the north and east, and with the Atlantic Ocean to the west and south. The North and Center regions of the Portuguese territory present a very mountainous terrain. The climate is predominantly temperate throughout the Portuguese mainland [20]. According to the Koppen classification, continental Portugal can be divided in two regions: one temperate climate with rainy Winter/dry and hot Summer, and the other with a temperate climate with rainy Winter/dry and not very hot Summer [21]. In terms of rainfall, its distribution is fairly asymmetrical throughout the territory, being more predominant in the North, moderate in the Center and scarce in the South, a region suffering by a progressive desertification process [22].

2.2. Sediment Delivery Ratio Model

The current soil erosion by water was modelled using the InVEST 3.6.0 software from Natural Capital Project [23]. InVEST models are “ready-to-use” spatially explicit models, that is, after the user collects and preprocesses the required input data, the model runs in a simple interface and delivers the expected outputs. The SDR model is based on the concept of hydrological connectivity requiring a minimal number of parameters [23]. The applied model uses the RUSLE (Revised Universal Soil Loss Equation) expression, where the factors are derived from different maps provided from different sources to determine the annual soil loss [23]. RUSLE is an extension of the original USLE (Universal Soil Loss Equation) with improvements in determining the factors controlling erosion [24,25]. This is an empirical model commonly used to estimate soil loss potential by water from hill-slopes across large areas of land. Soil-loss potential can be described as the estimated average annual value (ton ha⁻¹ y⁻¹) of soil erosion that an area can gather, using a factor-based approach with rainfall, soil erodibility, slope length, slope steepness and cover management and conservation practices as inputs [26]. Both USLE and RUSLE equations are written as follows [27]:

A = R∙K∙L∙S∙C∙P,

(1)

where A is the soil loss (ton ha^–1 y^–1); R is the rainfall erosivity factor (MJ mm ha^–1 h^–1 y^–1); K is the soil erodibility factor (ton ha h [ha MJ mm]^–1); L is the slope length factor; S is the slope steepness factor; C is the cover management factor; and P is the supporting practice factor; the L and S terms of the equation are often lumped together as “LS” and referred to as the topographic factor [27].

The software used to preprocess and analyse the geospatial data was ArcMap 10.7.1 for desktop [28]. All the input data had the ETRS_1989_TM06 coordinate reference system.

Table 1. Data sources for the data used as inputs for the Sediment Delivery Ratio (SDR) InVEST model.

Data	Source
Digital Elevation Model (DEM)	[29]
Rainfall Erosivity Index (R)	[30]
Soil Erodibility (K)	[31]
Land Use/Land Cover	[32,33]
P a and C b coefficients	[12,34]
Watersheds	[35]
Biophysical table	Created by authors of this study

^a Support practice factor, ^b Cover-management factor.

shows the data used as input for the SDR model in InVEST.

Relevant parameters used in SDR include the definition of the Threshold Flow Accumulation (TFA) values, which represent the number of upstream cells that must flow into a cell before it is considered part of a stream; two calibration parameters, k_b and IC₀, which determine the degree of connection from patches of land to the stream and percentage of soil loss that actually reaches the stream; and the SDR_max, which is the maximum SDR that a pixel can reach, in function of the soil texture. The default values were used, as indicated in the InVEST user guide documentation for this model [23].

The 30 m digital elevation model (DEM) was retrieved from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [29,36].

The rainfall erosivity index is an indicator of the ability of water to detach and transport soil particles; thus, erosion is sensitive to the intensity and duration of rainfall [26]. This index was provided by the Global Rainfall Erosivity Database (GloREDa) from the Joint Research Centre—European Soil Data Centre (ESDAC) and has a pixel resolution of 30 arc-seconds (~1 km at the Equador) [30]. GloREDa contains erosivity values estimated as R-factors from 3625 stations distributed in 63 countries worldwide. This is the result of an extensive data collection of high temporal resolution rainfall data from the maximum possible number of countries to have a representative sample across the different climatic and geographic gradient. It has three components: (i) the Rainfall Erosivity Database at European Scale (REDES) [37]; (ii) 1865 stations from 23 countries outside Europe; and (iii) 85 stations. Therefore, it is the most comprehensive global database including the largest possible number of stations with high temporal resolution rainfall data [38].

The soil erodibility factor (K-factor) is a lumped parameter that represents an integrated average annual value of the soil profile reaction to the processes of soil detachment and transport by raindrop impact and surface flow [24]. Consequently, K-factor is best obtained from direct measurements on natural plots [39]. However, this is a difficult task on a national or larger scale. To overcome this problem, measured K-factor values have been related to soil properties [39] estimating soil erodibility at the European level, based on attributes such as texture, organic matter, soil structure, and permeability, which were available from the Land Use/Cover Area frame Survey (LUCAS) [40] topsoil data [25]. Inverse distance weighting (IDW) was used to interpolate erodibility to a map with a grid-cell resolution of 10 km [6].

CORINE Land Cover (CLC) maps from European Environmental Agency (EEA) were the basis for the modelling approach [33]. CLC is a thematic land use/land cover cartography, available for the years 1990, 2000, 2006, 2012 and 2018, produced by the Directorate-General for the Territorial Development Portugal (DGT) for a project coordinated by the EEA. It consists of an inventory of land cover in 44 classes, with a Minimum Mapping Unit (MMU) of 25 hectares (ha) for areal phenomena and a minimum width of 100 m for linear phenomena [32]. The watersheds polygons were provided by the National Spatial Data Infrastructure (SNIG) [35].

The cover-management factor (C-factor) is used within both the USLE and the RUSLE to reflect the effect of cropping and management practices on erosion rates [34]. That is the most-used factor to compare the relative impacts of management options on conservation plans, indicating how the conservation plan will affect the average annual soil loss and how that potential soil loss will be distributed in time during construction activities, crop rotations or other management schemes [24]. The study made by Panagos and colleagues [34], in which the authors estimated C-factor values at a European level, was the starting point to estimate the C-factor values for the different land use/cover of the present study.

The support practices factor (P-factor) accounts for control practices that reduce the erosion potential of runoff by their influence on drainage patterns, runoff concentration, runoff velocity, and hydraulic forces exerted by the runoff on the soil surface. It is an expression of the overall effects of supporting conservation practices—such as contour farming, strip cropping, terracing, and subsurface drainage—on soil loss at a particular site, as those practices principally affect water erosion by modifying the flow pattern, grade or direction of surface runoff and by reducing the volume and rate of runoff [24]. The value of P-factor decreases by adopting these supporting conservation practices as they reduce runoff volume and velocity and encourage the deposition of sediment on the hill slope surface. The lower the P-factor value, the better the practice is for controlling soil erosion [12]. Many authors ignore the P-factor by giving it a value of 1, due to the difficulty of accurately mapping support practice factors [41]. In this study, the P-factor used for Portugal is 0.9178 for all CLC classes in the whole country, that is, the same value presented in the JRC-ESDAC study [12] for Portugal, under the assumption that all classes are arable lands [12]. Support practices over 0.95 have a greater influence in agricultural land, meaning that P > 0.95 will be prone to soil erosion.

The biophysical table (

Table 2. The biophysical table used in the SDR model, where ”LU-code” is the CORINE Land Cover (CLC) code for each land use class, “label” is the description of the class, and Universal Soil Loss Equation-c (“USLE-c”) and “USLE-p” are the cover-management (C)) and support practices (P) factors, respectively.

LU-Code	Label	USLE-c	USLE-p
111	Continuous urban fabric	0.1	0.9178
112	Discontinuous urban fabric	0.06	0.9178
121	Industrial or commercial units	1	0.9178
122	Road and rail networks and associated land	1	0.9178
123	Port areas	0.25	0.9178
124	Airports	0.25	0.9178
131	Mineral extraction sites	1	0.9178
132	Dump sites	0.9	0.9178
133	Construction sites	0.2	0.9178
141	Green urban areas	0.003	0.9178
142	Sport and leisure facilities	0.06	0.9178
211	Non-irrigated arable land	0.46	0.9178
212	Permanently irrigated land	0.36	0.9178
213	Rice fields	0.15	0.9178
221	Vineyards	0.4	0.9178
222	Fruit trees and berry plantations	0.3	0.9178
223	Olive groves	0.3	0.9178
231	Pastures	0.15	0.9178
241	Annual crops associated with permanent crops	0.35	0.9178
242	Complex cultivation patterns	0.2	0.9178
243	Land principally occupied by agriculture, with significant areas of natural vegetation	0.2	0.9178
244	Agro-forestry areas	0.13	0.9178
311	Broad-leaved forest	0.003	0.9178
312	Coniferous forest	0.003	0.9178
313	Mixed forest	0.003	0.9178
321	Natural grasslands	0.08	0.9178
322	Moors and heathland	0.1	0.9178
323	Sclerophyllous vegetation	0.1	0.9178
324	Transitional woodland-shrub	0.05	0.9178
331	Beaches, dunes, sands	0	0.9178
332	Bare rocks	0	0.9178
333	Sparsely vegetated areas	0.45	0.9178
334	Burnt areas	0.55	0.9178
411	Inland marshes	0	0.9178
421	Salt marshes	0	0.9178
422	Salines	0	0.9178
423	Intertidal flats	0	0.9178
511	Water courses	0	0.9178
512	Water bodies	0	0.9178
521	Coastal lagoons	0	0.9178
522	Estuaries	0	0.9178
523	Sea and ocean	0	0.9178

) was created using the CLC classes, and the C and P factors, as mentioned previously, by reviewing studies from the literature [12,34], and by adopting some values (e.g., for water bodies) from the biophysical table made available in the Natural Capital Project sample data [23]. In this table, the C-factor is represented by the USLE-c field, and the P-factor is represented by the USLE-p field. The LU-code field represents the CLC-code for each class.

The TFA values represent the number of upstream cells that must flow into a cell before it is considered part of a stream, which is used to classify streams from the DEM. Those are the values on which the model will create an output file of streams as close to reality as possible, depending on the DEM resolution, climate and topography [23]. IC₀ and k_b are two calibration parameters that determine the shape of the relationship between hydrologic connectivity and the sediment delivery ratio. The values for IC₀ and K_b are the default values for the SDR model, because it was determined that in hillslope-erosion specific yields, the model will perform better [42]. The SDR_max is the maximum SDR that a pixel can reach, which means that it is the maximum attainable sediment delivery ratio of the unit (in this case the pixel) [23,42]. The values for the SDR model are presented in

Table 3. Values used for the threshold flow accumulation, k_b, IC₀ and SDR_max parameters.

Parameters	Values
Threshold Flow Accumulation (TFA)	1000
kb	2
IC0	0.5
SDRmax	0.8

.

2.3. SDR Variation

The model output (Sediment Retention) with a spatial resolution of 30 m (equally as the DEM) was used for all the analysis. The expression used to calculate the sediment retention change between 1990 and 2018 was:

$$Sediment Retention Change \left(\%\right)=\frac{\left(S{R}_{2018}-S{R}_{1990}\right)}{S{R}_{1990}}\times 100,$$

(2)

where $S{R}_{1990}$ and $S{R}_{2018}$, are the raster outputs (Sediment Retention) from the SDR model, from 1990 and 2018, respectively.

2.4. Methodology

The sediment delivery retention model follows the workflow presented in Figure 2. Technical details of the InVEST SDR model can be obtained in the software user’s guide [23].

2.5. Model Validation

To validate the SDR model and its ability to assess soil erosion, output USLE was considered. This output represents the total potential soil loss by water per pixel in the original land cover calculated from the USLE equation [23]. A mean statistical test (t-test) was carried out to compare the mean results obtained for the NUTS III with our model and with the publicly available Soil Erosion by Water (RUSLE2015) dataset provided by European Soil Data Centre (ESDAC) [43]. The RUSLE2015 dataset uses a modified version of the RUSLE model, which delivers improved estimates based on higher resolution (100 m) peer-reviewed inputs of rainfall, soil, topography, land use and management from the year 2010 (i.e., the latest year for which most of the input factors are estimated) [12]. This dataset refers to the 28 Member States of the European Union, making it simple to extract the soil loss information for Portugal.

3. Results and Discussion

3.1. Main Results and Statitical Analysis

The SDR model was computed for five moments in time, corresponding to the years of the available CLC maps: 1990, 2000, 2006, 2012 and 2018. Along the 28 years evaluated, the sediment retention stayed fairly the same, ranging from 7.4 ton/ha in 1990 to 7.3 ton/ha in 2018, representing a decrease of 1.4%. The values for 2000, 2006 and 2012 were very similar, that is, 7.4, 7.3 and 7.4 ton/ha, respectively.

The SDR output maps for each of the years counted do not provide much information by themselves. Therefore, to better understand the outputs obtained, the raster calculator in ArcToolbox was used to calculate the percentage of gain/loss of sediment retention from 1990 to 2018. In Figure 3, it is possible to see that the difference of sediment retention throughout the territory is mainly between −5 and 5%, indicating that the territory did not suffer a big variation in terms of the capacity to retain sediments. Further analysis of the calculated raster shows the percentage of territory occupied by each class (

Table 4. Sediment retention change from 1990 to 2018 and the percentage of the total Portuguese territory occupation.

Class	Area per Class (km2)	Territory Occupation (%)
<−50	1314.95	1.21
−25–−15	1088.94	1.33
−15–−5	3972.48	4.85
−5–5	63,449.31	77.52
5–15	5557.27	6.79
15–25	2501.10	3.06
25–50	3726.13	4.55
>50	242.10	0.30
Total	81,852.28	100

). The results reveal that the sediment retention capacity is relatively the same throughout the Portuguese territory (77.52%) in the 28-year timeframe.

From 1990 to 2018, it is possible to observe that less than 1% of the territory had an increase in sediment-retention capacity and, oppositely, 4.85% registered a decrease. In Figure 3, this decline is mainly noticeable in the Alentejo regions (yellow colour). In those regions, the sediment retention capacity increased 4.55%. This may be possibly explained by seasonal erosion. A study concerning the Alqueva dam (located between Alentejo Central and Baixo Alentejo) explains that, according to the seasons, it is possible to verify different rates of soil erosion, in which the highest values occur during Autumn due to heavy rainfall, increasing the difference between rainfall erosivity and vegetation growth [44]. The GloREDa rainfall erosivity map used in this study has a high temporal resolution rainfall data (30–40 years), with mean values per year [38]. In this dataset, it is not possible to differentiate seasons that may result in discrepancies of sediment retention capacity in the same area. The authors of GloREDa rainfall erosivity map state that most of the uncertainty of this map is likely related to transition areas between different climatic zones and the different climatic conditions, ultimately resulting in high variability of rainfall amount, duration, magnitude and intensity [38]. A way to possibly overcome this issue is to create rainfall erosivity maps for each season for each year, and to use them in the InVEST SDR model, making the results more reliable. In this map, areas with a decrease of more than 50% in sediment retention capacity (i.e., 1.21% of the territory) can be found. When comparing these areas with the CLC for 2018 (Figure 1), it is possible to observe that many of them correspond to the Burnt area’s category, showing that burnt areas greatly influence the capacity of soil to retain sediments.

To understand which regions present a higher loss or gain in the capacity to retain sediments, a statistical analysis was applied to the map in Figure 3, using the zonal statistics tool from ArcGIS ArcToolbox. The map of Figure 4 shows the mean values differences (%) from 1990 to 2018 obtained per NUTS III after the classification in natural breaks. The regions represented in grey in the map of Figure 4 have fairly the same capacity of sediment retention throughout the years. Douro and the coastal regions are the ones that have a greater loss in sediment retention (peach colour), especially the region of Leiria (dark red colour), which was greatly affected by heavy forest fires in 2017. The Alentejo regions increased their capacity to retain sediments during the period of study (blue colour).

In Figure 5, it is possible to observe sediment retention (ton/ha) by NUTS III for each year. Alto Minho is the region with the highest capacity to retain sediments, while Lezíria do Tejo is the region with the lowest capacity. This figure also shows which regions have the highest variability in sediment retention per year. For instance, Cávado region has a high variability for the years 2000 and 2018; Beiras e Serra da Estrela region has high variability in 1990 and 2012; Aveiro and Médio Tejo regions both have a higher variation of mean values in 2012. These high levels of variability call for the need of further research to better understand their causes.

If wildfires directly influenced sediment retention losses, other causes that may justify the differences in sediment retention from 1990 to 2018 include changes in land use, especially for agriculture and urban growth. Another potential important explanation for the differences found in sediment retention is drought. According to the technical report of the European Environmental Agency [45], 2004/2005 was the year that suffered one of the worst droughts ever recorded in the Iberian Peninsula, with only half of the average precipitation, causing a decrease in the rivers’ flow. In 2003 and 2005, extreme fires followed by drought deeply affected the amount of sediment retention.

3.2. Model Validation

For the model validation, the model output USLE was used, the total potential soil loss by water per pixel in the original land cover calculated from the USLE equation [23]. A mean value was obtained for each of the 23 NUTS III regions for the year 2018 (

Table 5. Soil loss average value (ton/ha) for each NUTS III region in mainland Portugal, according to model output (USLE) for year 2018. Source: ESDAC dataset.

NUTS III	USLE	ESDAC (Reference)
Cávado	7.281	6.090
Ave	6.593	5.455
Área Metropolitana do Porto	4.351	4.455
Viseu Dão Lafões	3.593	3.256
Beira Baixa	2.186	0.980
Alto Tâmega	5.775	3.474
Tâmega e Sousa	8.742	7.643
Douro	11.859	6.039
Médio Tejo	1.996	0.866
Beiras e Serra da Estrela	4.165	2.761
Terras de Trás-os-Montes	4.910	2.716
Área Metropolitana de Lisboa	1.847	1.773
Alentejo Central	1.149	1.067
Algarve	2.206	1.871
Oeste	3.231	3.226
Região de Aveiro	1.476	1.320
Alto Minho	7.975	7.703
Alentejo Litoral	0.837	0.729
Baixo Alentejo	1.468	1.556
Região de Coimbra	3.689	1.312
Região de Leiria	3.984	1.013
Lezíria do Tejo	0.723	0.758
Alto Alentejo	1.305	1.052
Total (ton/ha)	67.117	91.340
Mean (ton/ha)	3.971	2.918

). Then, these values were compared with the ones using the ESDAC RUSLE2015 through a t-test. The null hypothesis was not rejected, that is, the observed difference of the sample means (3.971–2.918) was not enough to say that the means of USLE and RUSLE2015 differ significantly for the NUTS III regions. Thus, the model outputs are coherent with the ESDAC official data [12].

3.3. Limitations and Future Developments

The model produced some results that are estimates for the real world. In this sense, it must be taken into account that the results obtained will be partial and that the importance of the results obtained lies in the tendencies and insights that this analysis provides, which can be valuable in guiding the model development and making the best use of models of this type.

According to the InVEST models user’s guide, the SDR model presents some limitations [23]. The USLE [24] usage is very common, but this equation is limited in scope since it only represents rill/inter-rill erosion processes. Mass erosion processes such as landslides significantly impact to determine the amount of soil erosion in some areas. Nonetheless, those processes are not represented in this model. The SDR model is also very sensitive to k_b and IC₀ parameters, which are not physically based.

Another limitation is that the model produces NoData pixels in the stream network. The reason behind this is justified by the lack of in-stream processing. As it moves sediment down the slope, it stops calculations when the sediment reaches the stream, so in the estuary areas, where we have great water bodies, some pixel errors can occur in the water/land border. In addition, the SDR model is highly sensitive to most of the input data (due to its simplicity and the low number of parameters), which took a fair amount of time to process and adjust to the model. Additionally, the time it took to run and process the model, due to the heavy data inputs, was also a constraint.

Due to data availability limitations, climate-changing data was not considered. Additionally, for the rainfall data, the same JRC-ESDAC GloREDa map was used for all the temporal moments. This is an aspect that should be improved in future investigations.

Future developments should include a sensitivity analysis with advanced computational algorithms, such as Artificial Neural Networks, to determine how the model is affected when the values of the Borselli parameters k_b, the connectivity index IC₀, and the TFA values are calibrated to achieve the model’s optimal performance. Other future improvements should include the determination of the actual amount of sediments in each pixel to acknowledge where and how much soil gets deposited as it moves downhill towards a stream, or to quantify the erosion in the territory without converting the CLC classes as bare soil. A comparison of results from other types of models, such has the Soil and Water Assessment Tool (SWAT) [46] or the HEC-RAS [47], among others, should also be envisaged. Some studies have been made to determine the role played by USLE-based models on erosion rates (e.g., [41,48]), where different approaches to USLE-based models have been considered. Future development should be compared with the results of alternative modelling approaches. The use of different P-factor values can also be studied, doing an analysis per CLC class and applying it on the model for comparison purposes. Finally, further studies could also include the analysis of the main CLC classes per NUTS III and how these results can be correlated with soil erosion.

4. Conclusions

This study assessed the changes in sediment retention in mainland Portugal from 1990 to 2018. We quantified the effects of land use changes on the Portuguese hydrological basins and their impacts on soil erosion. Results show the different dynamics in sediment retention over the years at NUTS III level. The greater losses in sediment retention were observed in the Douro and coastal regions, and especially in the Region of Leiria. The model validation confirms that the outputs obtained are consistent with the ESDAC official data, demonstrating that the InVEST SDR model is an appropriate tool for estimating soil-loss potential by water at regional/national levels. Besides contributing with new information about sediment retention for Portugal in a 28-year timeframe, this study also provides a straightforward validation methodology of the results using credible reference datasets, which can be easily replicated for other study areas. The findings also contribute to the achievement of two Sustainable Development Goals (SDGs, Agenda 2030) [49]: Namely, goal number 15, “Life on land- Protect, restore and promote sustainable use of terrestrial ecosystems, sustainably manage forests, combat desertification, and halt and reverse land degradation and halt biodiversity loss”; and goal number 6, “Water and sanitation—Ensure availability and sustainable management of water and sanitation for all” [49].