Modelling sediment retention services and soil erosion changes in Portugal: A spatio-temporal approach

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 spatiotemporal 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 sediment retention capacity over the years. However, in most of the territory (77.52%) changes in sediment retention were little or not important (i.e., less than 5%). The statistical validation of the model proved the consistency of the results, highlighting the usefulness of this methodology to analyse the state of soil erosion in the country. These findings can be relevant to support strategies for more efficient land use planning regarding soil erosion mitigation practices.


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 are 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 partly through the detaching power of raindrops 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 to develop 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 worldwide [7][8][9][10][11]. At European level, [12] explored the use of the European Soil Erosion Model (EUROSEM) to simulate erosion processes, explicitly for rill and inter-rill flow. More recently, the RUSLE2015 model estimated soil loss at 100 m resolution for Europe [13]. A recent study conducted by [14] analysed soil loss and sediment exportation at the Winike watershed in Ethiopia, concluding that land use changes greatly affects the amount of soil loss in cultivated areas. Another recent study by [15] evaluated the soil erosion at a regional scale at Yunnan Province, China, using the Chinese Soil Loss Equation (CSLE) which allowed a more accurate soil erosion map for that province.
Particularly for Portugal, some studies have been carried out for modelling soil erosion at local scales (e.g. [16][17][18]). For example, [16] studied the nutrient retention by tradeoffs between sediments and vegetation types in Ria de Aveiro lagoon (central Portugal). [17] explored the effects of land abandonment on soil erosion and land degradation in the River Côa Valley (north-eastern Portugal). Recently, [18] investigated the influences of gully erosion in steep regions in the northern territory of Portugal. Albeit these studies have been made in different regions of Portugal, a deeper and validated study is yet to be carried to explain the effect of sediment retention on soil erosion in the entire territory. To fill this gap, the present study explores 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 between 1990 and 2018. Specifically, this study aims to: (i) estimate the soil loss at a pixel scale, and to (ii) estimate sediment retention variations at NUTS III level.
The study uses 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 with useful information to design a landscape effective planning for soil and water conservation.

Study Area
This study focuses on mainland Portugal ( Figure 1). Portugal, is a country in southern Europe, occupying a total area of 92,212 km2, whereas the mainland has a total area of 89,102.14 km2, with 23 statistical boundaries defined as NUTS III [19,20]. The mainland is located on the southwest of the Iberian Peninsula, bordering with 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 [21].

Sediment Delivery Ratio Model
The current soil erosion by water was modelled using InVEST 3.6.0 software from Natural Capital Project [22]. The InVEST models are "ready-to-use" spatially explicit models, i.e., after the user collect and pre-process 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 [22]. 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 [22]. RUSLE is an extension of the original USLE (Universal Soil Loss Equation) with improvements in determining the factors controlling erosion [23,24]. This is an empirical model commonly used to estimate soil loss potential by water from hillslopes across large areas of land. It estimates the annual soil loss that is due to erosion using a factor-based approach with rainfall, soil erodibility, slope length, slope steepness and cover management and conservation practices as inputs [25]. Both the USLE and the RUSLE equations are written as follows [26] (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 [26].
The software used to pre-process and analyse the geospatial data was ArcMap 10.7.1 for desktop [27]. All the input data had the ETRS_1989_TM06 coordinate reference system. Table 1 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, kb and IC0, 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 SDRmax, 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 for this model [22].
The 30 m digital elevation model (DEM) was retrieved from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [28,35].
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 [25] (Teng et al., 2016). This index was provided by the Global Rainfall Erosivity Database (GloREDa) from the Joint Research Centre -European Soil Data Centre (ESDAC) [29].
GloREDa contains erosivity values estimated as R-factors from 3,625 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) [36]; (ii) 1,865 stations from 23 countries outside Europe; and (iii) 85 stations collected from a literature review. Therefore, it is the most comprehensive global database including the largest possible number of stations with high temporal resolution rainfall data [37].
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 [23]. Consequently, K-factor is best obtained from direct measurements on natural plots [38]. 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. [38] estimated soil erodibility at European level, based on attributes (texture, organic carbon), which were available from the Land Use/Cover Area frame Survey (LUCAS) [39] topsoil data, using the original nomograph of [40]. Inverse distance weighting (IDW) was used to interpolate erodibility to a map with a grid-cell resolution of 10 km [6].
The land use/land cover products used in this project were the CORINE Land Cover (CLC) maps from European Environmental Agency (EEA) [32]. 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 [31]. The watersheds polygons were provided by the National Spatial Data Infrastructure (SNIG).
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 [33]. 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 [23]. The study made by [33], where 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 [23]. 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 [13]. According to [13], the P-factor used for Portugal is 0.9178 for all CLC classes.
The biophysical table (Table 2) was created using the CLC classes, and the C and P factors, as mentioned previously, by reviewing studies from the literature [13,33], and by adapting some values (e.g., for water bodies) from the biophysical table made available in the Natural Capital Project sample data [22]. 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. IC0 and kb are two calibration parameters that determine the shape of the relationship between hydrologic connectivity and the sediment delivery ratio. The SDRmax is the maximum SDR that a pixel can reach [22]. The values for the SDR model are presented in Table 3.

SDR variation
The model output (Sediment Retention) with a spatial resolution of 30m (equally as the DEM) was used for all the analysis. The expression used to calculate the sediment retention change between 1990 and 2018 was (2): where 1990 and 2018 , are the raster outputs (Sediment Retention) from the SDR model, from 1990 and 2018, respectively.

Model validation
To validate the SDR model and its ability to assess soil erosion, 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) (http://esdac.jrc.ec.europa.eu). 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) [13]. This dataset refers to the 28 Member States of the European Union, making it simple to extract the soil loss information for Portugal.

Results and discussion
The The SDR outputs for each of the years 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 between 1990 and 2018. In Figure 2, 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). The results reveal that the sediment retention capacity is relatively the same throughout the Portuguese territory (77.52%) in the 28 year's timeframe.  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 Fig. 2, using zonal statistics tool from ArcGIS ArcToolbox. The map of Figure 3 shows the mean values differences (%) between 1990 and 2018 obtained per NUTS III after the classification in natural breaks. The regions represented in grey in the map of Figure 3 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 the 2017 forest fires. The Alentejo regions increased their capacity to retain sediments during the period of study (blue colour).  In Figure 4 it is possible to observe sediment retention (ton/ha) by NUTS III for each year. Alto Minho is the region with better capacity to retain more sediments while Lezíria do Tejo is the region with the lowest capacity. 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 [41], 2004/05 was the year that has suffered one of the worst droughts ever recorded in the Iberian Peninsula, with only half of the average precipitation, causing the considerable decrease of the rivers flow. In 2003 and 2005, extreme fires followed by drought deeply affected the amount of sediment retention.

Model validation
For the model validation, the model output USLE was used. This output represents the total potential soil loss by water per pixel in the original land cover calculated from the USLE equation [22]. A mean value was obtained for each of the 23 NUTS III regions for the year 2018 (Table 5). Then, these values were compared with the ones using the ESDAC RUSLE2015 through a t-test. The null hypothesis was not rejected, i.e. 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 NUTSIII regions. Thus, the model outputs are coherent with the ESDAC official data [13].

Limitations
According to [22], the SDR model presents some limitations. The USLE [23] 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 kb and IC0 parameters, which are not physically based.
Another limitation is that the model produces NoData pixels in the stream network. The reason behind 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, it can occur some pixel errors in the water/land border. Besides, the SDR model is highly sensitivity 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 process the model, due to the heavy data inputs, was also a constrain.

Conclusions
This study assessed the changes in sediment retention in mainland Portugal between 1990 and 2018. We quantified the effects of land use changes on the Portuguese hydrological basins and its 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 frame, this study also provides a straight-forward validation methodology of the results using credible reference datasets. This methodology can be easily replicated for other study areas.
Future developments of this work should include a sensitivity analysis with advanced computational algorithms such as neural networks, to determine how the model is affected when the values of the Borselli parameters kb, the connectivity index IC0, and the TFA values are calibrated to achieve the model's optimal performance. Other future improvement 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 LULC classes as bare soil.

Data Availability Statement:
The data that support the findings of this study are available on reasonable request from the corresponding author.