1. Introduction
Soil is a critical natural resource that has been benefiting and supporting humankind in a variety of ways. Soil could be considered as a combination of parts of the geosphere, biosphere, atmosphere, and hydrosphere as it consists of solid mineral particles of debris from the weathering of rocks, organic material, and pores that contain air and water [
1]. Soil formation is affected by a few factors such as climate conditions, topography, the properties of the parent material, organisms (including both biota and human activity) and time. It is essential to know the rates of soil formation in watersheds to create sustainable watershed management strategies that will prevent soil degradation [
2]. It might take thousands of years just for a layer of a few centimeters of soil to develop [
3]. Time plays a vital role in the formation of soil as it takes many years for its generation but, on the contrary, soil loss is a process that happens considerably more rapidly. Special mention must be given to human activity, which influences the balance of the natural systems intensively through artificial interventions, leading to soil degradation, a worldwide concerning issue that deserves all of our attention.
According to the report «Caring for soil is caring for life» by the European Commission, soil degradation is quite evident in the EU as 60–70% of soils are currently maintained under unhealthy conditions, affecting not only agriculture but every aspect of our lives as the existence of humankind and life on Earth, in general, is strongly dependent on soil [
4]. Soil degradation occurs in many forms, such as soil erosion, pollution, salinisation, loss of organic matter and biodiversity, soil compaction, and surface sealing. Soil erosion is a phenomenon of great importance because of its intensity and variations in different areas. Especially in EU territories, it is estimated that approximately 1 billion tons per year of soil are lost because of erosion processes [
5]. For the reference year 2010, the mean soil loss rate in the EU was calculated using the RUSLE2015 model to be equal to 2.46 ton/ha/year which corresponds to 970 million tons of annual soil loss [
5]. Erosion is particularly prevalent in the Mediterranean basin due to extended periods of drought during summer followed by mild to wet winters with intense rainfall events which occur in areas with steep slopes and highly erodible soils, causing the loss of a significant amount of soil [
6,
7]. Greece is among the three Mediterranean countries (Italy, Spain, and Greece) showing the highest calculated soil erosion rates with 4.19 ton/ha/year which is higher than the corresponding value of 2010 which was equal to 4.13 ton/ha/year showing an increase of approximately 1.6% [
5,
8]. The establishment of soil erosion management practices is necessary for the maintenance of agricultural sustainability and the decrease of soil erosion rates. The control of soil erosion in areas of high risk can be achieved by implementing conservation practices such as the development and preservation of plant cover. In the study of Cedra et al., the biophysical impact of catch crops on soil erosion was investigated with rainfall simulation experiments showing that the catch crops caused a reduction of soil erosion in the examined plots [
9].
Reliable soil erosion data could be collected from field measurements. Despite the adequacy of this type of data, field erosion research is a high-priced process that would require a lot of time to create a database and ensure that the measurements are not affected by an extreme weather event or a few years of unusually high rainfall. Therefore, long-term measurements are required to study how soil erosion rates respond to land use and climate change. To study soil erosion, many methodologies have been developed by researchers, which intend to describe and quantify this natural phenomenon with the use of mathematical models, contemplating the factors that control the behavior of soil erosion through equations. Soil erosion modelling is the process of mathematically describing the detachment, transportation, and deposition of soil particles downstream [
10], which enables researchers to validate and codify the framework of erosion mechanisms [
11]. The fact that there are many empirical mathematical relationships for the estimation of soil erosion is because there is not yet a single relationship that could satisfy all possible erosion cases. In studies that focus on soil erosion monitoring in-site, the results from the analysis of the field measurements were compared with the calculated soil erosion rates generated from different models to validate and modify them [
12,
13].
The first attempt to develop an equation to calculate soil loss by correlating erosion rates with slope steepness and slope length on small fields belongs to Zingg [
14]. Some decades later, Morgan [
15] added specific factors for the improvement of this approach, such as the climatic [
16], the vegetation [
17], the conservation, and the soil erodibility factor. The change of the climatic factor in the rainfall erosivity factor eventually led to the Universal Soil Loss Equation (USLE) [
18]. The USLE and its adaptations, such as the Revised Universal Soil Loss Equation (RUSLE) [
19] and the Modified Universal Soil Loss Equation (MUSLE) [
20], are the most frequently utilized soil erosion models which calculate the average annual soil loss due to sheet and rill erosion. Moreover, for the European territories, the European Soil Erosion Model (EUROSEM) [
21] and the Pan European Soil Erosion Risk Assessment (PESERA) [
22] were generated.
The Revised Universal Soil Loss Equation (RUSLE) is the most widely used model that calculates the rate of soil erosion induced by rainfall and overland flow. The RUSLE model is widely used in Greek territory by researchers, with the necessary adaptations in relation to the various conditions of each research area [
5,
7,
23,
24,
25,
26,
27,
28,
29]. Its ease of use and the ability to define the variables causing water erosion while using highly accessible and accurate data to calculate the average annual soil loss rates in each study area are the elements that have contributed to its widespread adaptation. The factors that contribute to the RUSLE model are the Rainfall-Runoff Erosivity factor (R), the Soil Erodibility factor (K), Slope Length and Steepness factor (LS), Cover Management factor (C), and the Support Practice factor (P).
The accuracy of the input data and the methodologies used to compute each factor have a direct impact on the quality of the final computation. Specifically, the LS factor is of great interest because of the variety of applied equations in the literature to calculate it and the fact that these equations require the use of Digital Elevation Models (DEM). Τhe use of DEM for describing the topography of each research area is a parameter that has a substantial impact on the soil erosion modelling outcomes as the DEM source and resolution determine the accuracy of the derived topographical characteristics like elevation, slope, and slope length [
30]. Concerning the definition of the Digital Elevation Model, it is a broad term for describing the digital and mathematical representation of topographic surfaces [
31]. In the case where the digital model is referring to the bare ground surface, it is called a Digital Terrain Model (DTM), but if the surface includes any feature above the ground, such as trees and buildings, then it is called a Digital Surface Model (DSM) [
32,
33].
Many researchers have applied different methodologies for the calculation of the LS factor as well as different types of elevation data like DTMs and DSMs. In 2015, Panagos et al. acquired a 25 m resolution DEM for the estimation of the LS factor for the European Union, indicating that DEM resolution does have a significant impact on the spatial distribution of the LS factor since the increase of DEM resolution contributes to the more accurate representation of the landscape and leads to the improvement of the spatial soil loss calculations [
34]. Efthimiou et al., examined the effect of topography on soil erosion by calculating the LS factor using the proposed equations from the USLE and RUSLE models in 8 sub-watersheds in northwestern Greece and pointed out that estimated soil erosion from the two models matches for slopes <20%, but for steeper slopes above 20% soil loss from USLE calculated higher values than RUSLE [
35]. Fu et al., used DEM with different grid sizes from 2 to 30 m to calculate the LS factor with the Chinese Soil Loss Equation (CSLE) and concluded that both the average LS factor and soil loss decreased with the increase of the DEM grid size, although the corresponding values from the 10 m and 2 m were similar [
36]. Mondal et al., performed an analysis on open source DEM with 30 m resolution and aggregated grid sizes (90, 150, 210, 270, and 330 m) regarding their vertical accuracy and uncertainty in soil erosion estimation using the RUSLE model, indicating that both the accuracy of the elevation data and the estimated soil erosion decrease when the DEM grid size increases [
37]. In contrast to the previous studies, Shan et al. used DEM with resolutions from 1 to 90 m and the RUSLE model to examine the effect of resolution on the computations of LS factor and soil erosion in a mountainous area in Australia and suggested that there is an increase in the LS values with the increase of the DEM grid sizes and, more specifically, the L factor increases and the S factor decreases [
38]. These opposite results may be linked to the methodology of the LS factor, which did not consider a cutoff factor of the slope length, which would cause the overestimation of the LS factor in areas with long slopes where runoff concentration would produce gullies [
39] and to the morphology of the study areas, which were characterized by larger relief regions compared to other studies [
40]. In the study of Kruk et al., high-resolution DEM were generated from different methods (aerial photographs, aerial laser scanner, and terrestrial laser scanner) to assess the LS factor and soil erosion using the USLE model, and it was shown that resolution has a significant impact on the LS factor, with the S coefficient strongly affecting the calculation of LS [
41]. Lu et al., derived topographic data from five different sources with spatial resolution from 5 to 90 m to assess the LS factor in five catchments in China, and the results indicated that by increasing the grid size, the S factor decreases, the L factor increases, and the overall accuracy of the estimated LS factor is decreased [
39]. Wang et al., studied the effect of DEM resolution on the LS factor by implementing two DEM, one derived from Lidar datasets with a 5 m grid size and the other obtained from an SRTM dataset with a 30 m grid size for the computation of LS in 24 watersheds in North America [
40]. This study highlighted the increase in the L factor and the decrease in the S factor with coarser DEM resolutions and, regarding the LS factor, the results showed an increase in areas with large relief and a decrease in smaller relief areas. Azizian and Koohi focused on the impact of different source DEM (ALOS, ASTER, SRTM) with different resolutions from 30 to 200 m and produced LS factors with various equations as well as the sediment yield rates in a study area in Iran, concluding that DEM resolution has a more vital role than DEM source and the different methods for the calculation of LS factor [
42]. Another study that assessed the impact of DEM from different sources (ASTER, SRTM, and CARTOSAT) and various resolutions from 30 to 250 m was carried out by Pandey et al., who performed the RUSLE model for the Karso watershed in India and displayed the scale dependence of RUSLE [
43]. Regarding the effect of DEM resolution, in the above study there is an increase in the calculated mean and maximum LS factor with an increase of the grid size up to 50 and 100 m but for higher grid sizes a corresponding decrease of these values is noticed. Moreover, the impact of the DEM source highlighted the weaker performance of ASTER DEM compared with SRTM and CARTOSAT. Kumar et al., applied the RUSLE model and acquired DEM from SRTM, ALOS, and MERIT with a resolution of 30 and 90 m and compared the calculated soil erosion with observed data, which led to the conclusion that the DEM with a finer resolution (30 m) performed better than the coarser ones (90 m) [
44]. Fijałkowska executed a comparative analysis on elevation, slope, L factor, S factor, and total LS factor values using DTMs with various resolutions (1 m to 90 m) and from different sources using the RUSLE model for 3 study areas in Poland, emphasizing the non-uniformity of the change of LS factor with the increase of DTM overall accuracy and concluded that high resolution and accuracy data should be assessed in soil erosion modelling [
45].
The current study focuses on the importance of topography in soil erosion modelling by examining the impact of topographic data from various sources in the calculation of the slope length and slope steepness factor (LS) of the RUSLE model. It is particularly applied in an area with dense morphology, which has suffered extended surface changes due to repeated wildfires, and this is the drainage basin of the Pinios earth-filled dam located in the Ilia Regional Unit, Western Greece. For the purposes of this study, the topographic data included elevation data from 4 different sources that were compared with Ground Control Points (GCP) to assess their relative vertical accuracy and suitability to produce accurate soil erosion models. Specifically, six Digital Elevation Models with various resolutions and two different methods were applied for the calculation of the LS factor in the examined area. Moreover, for the estimation of the LS factor, the equations proposed by Mitasova et al. [
46,
47] and the CALSITE (Calibrated Simulation of Transport Erosion) model [
48,
49] were used.
This research aims to investigate and derive the correlation between DEM and the LS factor, as it is the most crucial parameter of the RUSLE model for estimating soil erosion annual rates. The inspection of the quality of the used DEM in the examined area is of particular interest due to its dense morphology and susceptibility to soil erosion, strongly related to the major wildfires that have affected it for the last fifteen years.
The methodology applied in this research consists of four main stages: (a) the comparison of the examined DEM elevation values with the corresponding values derived from GCP; (b) the comparison of the DEM slope values; (c) the calculation and comparison of the LS factor using two methods, which are referred to as the MIT method [
46] and the CAL method [
48,
49], and (d) the estimation and comparison of the annual soil loss in the pilot area using the RUSLE model and the different calculated LS factors.
5. Discussion
The evaluation of role of topography in soil erosion modelling was assessed by examining the impact of elevation data from various sources on the slope length and slope steepness factor (LS) estimation. Specifically, 6 Digital Elevation Models (DEM) were acquired from 4 different sources and then compared with ground control points (GCP) to analyze their relative vertical accuracy and evaluate them in terms of their suitability for soil erosion model production. The DEM had various resolutions and were used to calculate the LS factor with two different equations. The RUSLE model was then used to estimate the annual soil erosion rates for each one of the examined approaches and DEM. Regarding the soil loss estimation with the RUSLE model, it is noted that all the parameters (R, K, C, and P) remained constant, except for the LS factor, which was calculated by applying two different equations, meaning that any variations in estimated soil loss values for each DEM are caused by the LS factor.
The quality of the final computations of the LS factor, and subsequently the annual soil erosion rates, directly depends on the accuracy of the input data and the calculation techniques utilized for each factor of the RUSLE model. In particular, the LS factor is very important because of the variety of equations used in literature for its estimation and the fact that it requires the application of DEM data. The accuracy of the DEM-generated topographical characteristics, such as elevation, slope, and slope length, is determined by the DEM source and resolution, which has a significant impact on the soil erosion modelling results as the DEM are used to describe the topography of each region of interest.
According to Polidori and Hage, there are a few quality criteria to examine the adequacy of a DEM for a specific use, which could be divided into two sectors, elevation and shape/topologic quality [
32]. The first type refers to the absolute or relative accuracy of a DEM, and the second one corresponds to the deliverables of a DEM which describe the relief, such as slope and aspect. DEM quality assessment refers to both vertical and horizontal accuracy, even though in many case studies only the vertical component is analyzed with the calculation of statistical indices.
In this research, 6 DEM (2 DTM and 4 DSM) from different sources and with various resolutions, were acquired to evaluate their suitability for soil erosion modelling. The first section of this study contains the evaluation of the relative vertical accuracy of the examined DEM by comparing them to GCP. Initially, the error (Δh) was calculated by subtracting the altimetric value of each checkpoint from the corresponding value of each DEM. The basic statistics of the error, as well as more sufficient indices like MAE, MSE, RMSE, and sMAPE, were calculated too. These metrics showed that the DTM derived from the Hellenic Cadastre (DEM5) and the corrected one (DEM5_IDW) are performing better than the other DEM with RMSE values of 1.80 m and 0.54, respectively. Following this, ALOS30 and SRTM30 appear to have less satisfactory RMSE values equal to 8.95 m and 11.67 m, respectively, while ASTER30 (RMSE: 18.27 m) and SRTM90 (RMSE: 20.02 m) are underperforming with RMSE values of approximately 20 m. The calculated RMSE value of ALOS30 for the pilot area is also supported by the studies of Nikolakopoulos [
33] and Stamatiou et al. [
65] which reported RMSE values for the Greek territory equal to 8.58 m and 8.75 m, respectively. The RMSE value of ASTER30 was quite higher, almost double, the corresponding value that was estimated by other studies for the Greek territory, which reported an RMSE value equal to 10 m [
68,
69]. That is also the case for SRTM30, as reported RMSE values from other studies equal to 6.1 m [
69] and 6.8 m [
69] are almost half of the calculated one for the pilot area. The thorough eminence of the DEM5 and DEM5_IDW is also endorsed by the estimation of two more robust metrics, the MAD and NMAD. These metrics follow almost the same sequence from lowest to highest value as RMSE, with the exception that ASTER30 appears to have the highest MAD and NMAD values, followed by SRTM90. The underperformance of ASTER30 may be linked to the fact that the study area is mainly occupied by agriculture, which occurs in areas with relatively flat terrain where ASTER DEM shows lower vertical accuracy [
68,
73].
The comparison of DEM with different cell sizes provides essential information regarding the accuracy of their elevation data, while the comparison of the DEM derivative files such as slope is a more complicated process. The similarity of the statistics of the DEM elevation values is quite high even though they have different resolutions and sources, but in the case of a slope, the effect of the resolution is more evident. According to the work of Carrera-Hernandez, who examined the vertical accuracy of 8 DSMs from different sources in Mexico, it was shown that slope is the major factor that influences the error of DSMs among vegetation and aspect [
78]. Zhou et al. analyzed the errors of slope and aspect with the used algorithms and the DEM properties and pointed out that even though DEM with high resolution have more details, the data error tends to increase, in contrast to the low-resolution DEM in which the data error has less influence while the error associated with algorithms increases, meaning that high-resolution DEM does not guarantee high accuracy for its deliverables [
108]. The comparison of slopes for the evaluated DEM in the research area showed considerable variations. The DEM resolution has a major impact on the calculation of slope. The frequency of steeper slopes tends to increase in DEM with high resolution, as the DEM5 and DEM5_IDW with a 5 m resolution have the highest maximum and mean values, compared to the SRTM90 with a 90 m resolution, which shows the lowest maximum and mean values of all DEM.
DEM resolution and DEM source are very important parameters that influence the accuracy of soil erosion models [
42,
43]. Most of the time, the DEM have been used in soil erosion modelling studies without considering the examination of their adequacy regarding the accuracy and the quality of the source data which plays a significant role in the outcomes [
30,
45]. It is also important to consider the fact that both DTM and DSM data have been used in soil erosion modelling. By definition, a DTM refers to the bare ground surface, and it is considered to be the correct type of DEM for applications such as soil erosion simulations. The implementation of DSM data in soil erosion modelling is mainly due to the wide availability of freely acquired DSMs like those that were used in the current research. The fact that DSMs include the elevation of the canopy and man-made constructions can produce substantial errors in the calculations of the LS factor and, subsequently, soil erosion.
The effect of DEM in soil erosion modelling was examined by many researchers using the LISEM model (Limburg Soil Erosion Model) and DEM with different pixel sizes and concluded that large pixel sizes tend to cause a reduction in the estimations of soil loss, which is influenced primarily by the reduction of the DEM derived slope caused by the usage of DEM with large pixel size [
42,
113,
114]. In the study of Chidi et al., the RUSLE model was applied using DEM with varying resolutions from 5 cm to 10 m in the Middle Hill area in Nepal [
115]. This research pointed out that by increasing the pixel size in DEM, the expected slope gradient will decrease, influencing the LS factor and the estimated soil loss, which also decreases with the increase of DEM pixel size [
115]. The current statement that LS factor values decrease with coarser DEM is following the studies of [
36,
37,
39,
115] and is in contrast to the study of Shan et al. [
38].
The high LS values of the MIT method are related to the fact that this LS equation is an integral part of the USPED model (Unit Stream Power—based Erosion Deposition), and therefore should be interpreted in the context of this model. It should also be considered that the USPED model [
47] allows the estimation not only of the net erosion but also determines the deposition in a basin, unlike the RUSLE model [
116]. According to Karásek et al., the LS factor calculated with this method is characterized by a line structure of pixel cells in which the LS factor values are increasing accordingly to the direction of overland flow downstream of the slope [
117]. This leads to a sharp transition between the high and low LS values, which are associated with the use of the flow accumulation in this equation and result in quite high mean LS values. An important difference between these approaches is that the MIT method uses the upstream contribution area or the specific catchment area (A) whereas the CAL method uses the slope length (λ) for the calculation of the slope length or L factor. More specifically, in the case of the CAL method, a threshold of 305 m was considered in the slope length calculations. Regarding the slope steepness, or S factor, a threshold slope value of 9% was taken into consideration. It is also noted that the calculated mean LS factors from the MIT and the CAL method were compared with the estimated mean LS value that was calculated by the European Soil Data Centre (ESDAC), which developed a pan-European high-resolution soil erosion assessment. This comparison was conducted in the circumstance that there are no field measurements available for the evaluation of the LS factor. This LS factor assessment showed that the LS values calculated from the CAL method, with a range of 3.24 (from DEM5_IDW) to 5.08 (from SRT90), were very close to the corresponding value from ESDAC, which is equal to 3.31, while the corresponding values from the MIT method are significantly higher, ranging from 8.77 (from SRTM30) to 11.41 (from ALOS30). This indicates that the CAL method performs better in the pilot area.
6. Conclusions
The objective of this research is to examine the connection between DEM and the LS factor, which is one of the most significant parameters of the RUSLE model, to estimate soil erosion annual rates. Due to the dense morphology of the research area and its susceptibility to soil erosion, the examination of the quality of DEM has proved to be an important issue.
Regarding the DEM analysis, which was performed in this study, the elevation data showed a satisfactory similarity based on the high values of the calculated Pearson correlation index. Based on the statistical assessment, the DEM5 and the corrected DEM5_IDW were found to be better compared with the other DEM, followed by ALOS30, SRTM30, ASTER30, and SRTM90, and the calculated error increased with the increase of the elevation (>1000 m) for all DEM. The comparison of the slopes from each of the 6 examined DEM showed a greater discrepancy, mainly because of the different DEM cell sizes and their types, as both DTMs and DSMs were acquired for this research.
The computation of the LS factor was performed with the application of two different methods. The first method proposed by Mitasova et al. [
46,
47] showed much higher LS values than the second one proposed by the CALSITE method. Increasing the DEM grid size, the calculated LS mean values were increased in the case of the CAL method, but for the MIT method, the results did not follow a specific sequence. Regarding max LS values, the DEM with the 5 m resolution showed the highest value in both cases. The lowest value was calculated from the SRTM90 with the MIT method, in contrast to the CAL method, as the SRTM90 had the second highest value. The evaluation of the correlation using the Pearson coefficient showed that the LS values calculated with the CAL method for each DEM have a higher similarity with each other than the LS values for each DEM from the MIT method, which showed a great divergence. Both of the LS factor computation methods presented significantly higher values for slopes that exceed 20° for all the examined DEM.
As for soil loss, the mean SE annual rates calculated for each DEM using the CAL method are significantly lower than the relevant SE values calculated using the MIT method, following the same pattern as the calculated LS factor. The Pearson correlation index showed greater similarity in SE values calculated with the LS factor from the CAL method compared to the MIT method, and in both cases, SE values estimated for slopes above 20° were substantially increased.
The different results derived from various studies regarding the LS factor are related to the different approaches in each applied methodology, the special characteristics related to the landscape for the area of interest, the different sources of examined DEM, and the different LS factor equations. The comparison of the LS values calculated with the two examined approaches and with the use of different DEM data with various resolutions and from diverse sources does not change consistently with the increase of DEM grid size and accuracy. This is also supported by the study of Fijałkowska, which stated that the relationship between DEM resolution and LS factor is indeed complex because the LS factor from one DEM can produce high values for the same area where the LS factor from another DEM calculates low values [
45].
The use of adequate data is essential for soil erosion modelling. In particular, the accuracy of topographic data like DEM plays a vital role in the identification of soil erosion high-risk areas. The precise detection of soil erosion-prone areas is very important for protection measurement planning. The sustainability of the examined areas defines the types of applications that need to be fulfilled. Nature-based solutions are considered to be proper, mainly for areas occupied by a high percentage of agricultural land, as their benefits are long-term and cost-efficient [
118].
In conclusion, this research highlighted the importance of evaluating the primary topographic data and the different approaches to LS factor calculation before continuing with the assessment of soil erosion modelling, especially in areas with diverse landscape. It seems that the use of an LS equation that imports thresholds in its formula, such as the CALSITE method, helps to avoid overestimation in soil loss calculations. The next step requires field measurements to evaluate and validate the different data and methods used regarding the LS factor and to define which method is the most appropriate for the corresponding area of interest. Towards this aspect, different methodologies and models like WATEM/SEDEM and the InVEST SDR model could be acquired.