Analysis of the Influence of DTM Source Data on the LS Factors of the Soil Water Erosion Model Values with the Use of GIS Technology

: Counteracting soil degradation is one of the strategic priorities for sustainable development. One of the most important current challenges is effective management of available resources. Multiple studies in various aspects of soil water erosion are conducted in many research institutions in the world. They concern, among others, the development of risk estimation models and the use of new data for modelling. The aim of the presented research was a discussion on the impact of the accuracy and detail of elevation data sources on the results of soil water erosion topographic factors modelling. Elevation data for this research were chosen to reflect various technologies of data acquisition, differences in the accuracy and detail of field forms mapping and, consequently, the spatial resolution of the digital terrain models (DTMs). The methodology of the universal soil loss equation USLE/RUSLE was used for the L and S factors modelling and calculation. The research was carried out in three study areas located in different types of geographical regions in Poland: uplands, highlands and lake districts. The proposed methodology consisted of conducting detailed comparative elevation and slope value assessments, calculating the values of topographical factors of the universal soil loss equation: slope length (L) and slope (S) and a detailed analysis of the total LS factors values. An approach to assess LS factors values within homogeneous areas such as agricultural plots has also been proposed. The studies draw the conclusion that the values of topographical factors obtained from various DTM sources were significantly different. It was shown that the choice of the right modelling data has a significant impact on the L and S factors values and, thus, also, on the decision-making process. The conducted research has definitely shown that data of the highest accuracy and detail should be used to study local phenomena (inter alia erosion), even analysing a large area.


Introduction
Soils are the basis for plant production and food security. The 2015 Food and Agriculture Organization (FAO) report "Status of the World's Soil Resources" began with the statement: "Soils are fundamental to life on Earth" [1]. In paragraph 10, The Updated FAO World Soil Charter of 2015 emphasised that soil degradation significantly reduces or even eliminates the natural functions of soil. Preventing or minimising soil degradation is necessary to maintain basic functions of soils and is much more economically viable than the reclamation of degraded soils [2,3]. Unfortunately, it is estimated that most of the world's soil resources are in merely adequate, poor or very poor conditions, and 33% of the land is moderately-to-severely degraded due to, among others, water, wind and crop erosion. It is predicted that a further loss of soil productivity threatens the balance and security of food production, which will cause further price increases and potentially worsen the quality of life of millions of people [2].
In 1983, Renard and Foster, working on a mathematical description of the factors causing water erosion, proposed a general formula showing the dependence of the extent of soil loss on the factors that cause it [4]. Since then, various approaches to erosion modelling have been developed (expert-based [5,6], conceptual, empirical [7][8][9] and physicsbased/mechanistic models [7,10]), and these different models have different versions that have undergone many changes and improvements as research has developed [11]. Empirical models constitute an important group of erosion models. These formulas are developed on the basis of the results of research conducted on experimental plots. It is estimated that, for calibrating empirical models of water erosion, we currently have results from over 24,000 measurements conducted mainly in the USA, Europe, Brazil and China [12]. Unfortunately, they are often impossible to compare due to the large diversity of climatic zones and the characteristics of individual slopes [7][8][9]13,14]. Measurements on experimental plots of a given size and with an area limited by natural terrain barriers are not able to adequately reflect the dynamics of the erosion process. During each erosion event, the distance of displacement of the eroded soil material may be different [14]. Moreover, this type of research requires significant financial outlays and is time-consuming, and the obtained results often cannot be extrapolated to larger areas [7,9,15,16]. Despite the above-mentioned limitations, the ease of implementation of empirical equations and the availability of data make these models one of the most popular methods used to assess the impact of water erosion.
The most widely used model in the empirical model range is the Universal Soil Loss Equation, developed in 1965 [17] and its modifications presented in many later studies. The formula for calculating the average annual soil loss is presented in Equation (1): A= R * K * C * L * S * P (1) where: A = mean annual soil loss (in t ha −1 yr −1 ), R = rainfall erosivity factor (in MJ mm ha −1 h −1 yr −1 ), K = soil erodibility factor (in t h MJ −1 mm −1 ), C = cover management factor (dimensionless), L = slope length factor (dimensionless), S = slope steepness factor (dimensionless) and P = support practice factor (dimensionless). Models based on the universal soil loss equation USLE/RUSLE methodology [18] are commonly used to estimate the risk of soil water erosion in analyses covering small areas (single farmland and slopes), as well as in analyses for large areas (one or more catchments). Some of the research conducted was undertaken by . The European Soil Data Centre (ESDAC) and the European Soil Bureau (JRC) also use these equations when calculating soil loss, including in the PESERA, G2 and RUSLE2015 models [42]. Various authors have emphasised the importance of proper determination of the impact of each factor included in a given erosion model: the precipitation [25,49], the natural susceptibility of soils to erosion [22,50] and land cover and land use [42,50,51], but the vast majority indicate that the factors of length (L) and slope (S) are those whose precise calculations determine the obtained results to the greatest extent. Additionally, uncertainty analyses indicate the significance of these factors [19,22,29,36,37,42,[49][50][51][52][53][54][55][56][57].
The models mentioned above differ, inter alia, in the formulas that are used to calculate the impacts of individual factors. In determining the impact of the topography, the most common calculation of the L factor value is based on the distance travelled by the runoff water or on the basis of the unit upslope contributing area, understood as the area from which the water flows to a given point [19] (considered to better reflect the erosive impact of the track passed by water flowing down the slope than a simple linear flow path). The determination of unit upslope contributing areas involves the identification of flow directions for this purpose and the calculation of flow accumulation on its basis.
Algorithms that can be used to determine the paths of water flowing down the slope is the assumption of water runoff from the higher raster cells to only one cell (single flow) or to several lower cells (multiple flow). The most frequently used single-flow algorithms are D8, Rho8 and KRA, while the multiple flow algorithms include, among others, DF8, FRho8, D∞, MD∞ and DEMON [37,58]. Despite the expected accuracy of the obtained results, the tests of the above-mentioned algorithms do not always yield similar results, which is often caused by the complexity of the terrain being analysed [59]. Some authors emphasise a much better efficiency of the DEMON algorithm or divergent runoff models [23,37,60], while others indicate no significant differences in the course of the hydrographic network determination and the value of the runoff accumulation [61,62]. The common use of the USLE/RUSLE models facilitates a comparison of the obtained results. A review of the literature sources showed that various data are used to calculate water erosion soil losses and to indicate areas particularly exposed to it (so-called "hot spots"). It can be seen that the obtained results were sometimes significantly different. This is confirmed by the comparison (review) of the selected results of European studies presented in Table 1.
The analysis of the above list led to the following questions: why are there such differences in the obtained results (maximum soil loss values vary from 6.4 t ha −1 y −1 to 11,680 t ha −1 y −1 , while research at the pan-European level reaches little more than 50 t ha −1 y −1 )? Where do they come from? What causes such discrepancies? Which factors contribute to such differences? If it is the influence of the length and gradient of the slope, do the obtained differences result from the digital terrain model (DTM) source data, their resolution, the formula used or its implementation in a software used? Table 1. Review of soil loss results from selected European studies in the period 1996-2020, and a comparison with the results of field studies from USA and Europe investigations. Equation: U (Universal Soil Loss Equation (USLE)) and R (Revised Universal Soil Loss Equation), flow: flow length (l) and unit contribution area (a) and flow algorithm: single flow direction (SFD) and multiple flow direction (MFD). The literature review gives the impression that, usually, a widely available data is used without first considering the appropriate accuracy and detail of the source data and their impact on the modelling results. Very often, a medium-scale DTM is used to model a phenomenon that is clearly local in nature. This approach is understandable for analyses performed in the absence of available data of adequate detail and accuracy, but it raises doubts in times of the widespread availability of high-resolution elevation data. Although the term "high-resolution DTM" has significantly changed its meaning since the inception of Geographical Information System (GIS) technology, many authors dealing with soil water erosion indicate that the use of source data with high accuracy and detail is an essential aspect of erosion modelling [34,37,55,[83][84][85]. Moreover, some studies concerning the uncertainty of the model results indicate that the uncertainty of estimating the values of the L and S factors is crucial. This impact can be reduced by using DTM with very high accuracy and detail [11,39,52,[85][86][87][88]. Initially, the data sources for calculating the values of the L and S factors came from field studies; then, they were obtained by vectorisation of the base of topographic maps in different scales [39] and, then, from DTM from various sources (active and passive satellite data [38,56], aerial photography or laser scanning [85,89,90] [42], as well as models from aerial photos and ALS [85,89,90]. For individual fields, GNSS measurements were also made in a regular grid [60]. In the project implemented for the entire EU area, the values of the L and S factors were calculated from the EU-DEM model (SRTM + ASTER) with a spatial resolution of 25 m, as it is the highest spatial resolution of a homogeneous DTM available for all of Europe. Individual Member States have used DTM with a resolution of 10-40 m to model the risk of water erosion [93].

No
However, the sensitivity of the erosion models to the source data used has not been comprehensively investigated so far, although the publications emphasise the impact of the accuracy and detail of the source DTM on the values of the L and S factors obtained as a result of its processing [39,43,94]. Tests evaluating the effect of reducing the spatial resolution of source DTM on the modelling results were conducted by, among others, Molnar and Julien [95], Wilson [58], Lee and Lee [88], Liu et al. [96] and Zhao et al. [83].
The methodology used by these researchers assumed the study of the dependence of the obtained results on the spatial resolution of DTM and the subsequent versions of the original DTM, with increasingly lower resolutions (obtained as a result of degradation of the source data resolution). Such models are consistent and only simulate the real data with lower resolutions. However, such comparisons did not make possible the check of the impact of the real data source on the erosion modelling (e.g., possible shifts of terrain forms, lack of visibility of minor terrain forms; this problem is well-reflected in comparative analyses of the hydrographic network system and catchment boundaries obtained from DTM from various sources). Taking into account the high availability of various elevation data for the area of the country or even on a continental scale (obtained with the use of different technologies and characterised by different accuracies and details). There is a need to conduct in-depth research aimed at determining the impact of the source of elevation data on the results of determining the values of the L and S factors of soil water erosion modelling and, consequently, on the determination of areas at risk of water erosion.

Study Areas
The proposed methodology was applied in 3 research areas, the selection of which was preceded by an analysis of data available for Poland describing the risk of water erosion developed by The Institute of Soil Science and Plant Cultivation (IUNG). The water erosion risk database distinguishes 6 levels (an expert classification based on: slope class, soil texture classes and total annual rainfall class) of erosion risk: from 0-no risk to 5very strong erosion risk (dimensionless). The research areas are characterised by varied topography and were selected in such a way that they contain areas of high erosion risk: levels 4 and 5, as well as areas with no risk and low erosion risk: levels 0 and 1. The research areas are located in various types of physical and geographical regions of Poland: foothills, highlands and lake districts ( Figure 1). The research area of Janowice covers an area of 1355 ha and is located in the Western part of Pogórze Ciężkowickie. The terrain heights (based on ALS data) range from 199 m above sea level up to 454 m above sea level. The maximum slopes of the terrain are significant and reach values up to 78° (ALS data) ( Figure 2a).The land is mainly covered by grasslands and agricultural crops (55.9%), forests and wooded areas (29.3%), built-up areas (5.8%) and permanent crops (4.5%) (based on the BDOT10k topo reference database). The research area of Winiarki covers an area of 277 ha and is located in the eastern part of the Sandomierz Upland. Terrain heights (ALS data) range from 138 m above sea level up to 208 m above sea level, and the slope values (ALS data) in this area reach 74° ( Figure 3a). The land is mainly covered by permanent crops (38.1%), grasslands and agricultural crops (32.5%), forests and wooded areas (23.3%) and built-up areas (5.5%) (based on the BDOT10k topo reference database). The research area of Bolcie covers an area of 501 ha and is located in the northern part of the East Suwałki Lake District. The north-western part of the area border coincides with the Polish and Lithuanian borders. The terrain altitudes (ALS data) vary from 206 m above sea level up to 273 m above sea level. The slope values in this area (ALS data) reach up to 59° (Figure 4a). The land is mainly covered by grasslands and agricultural crops (89.9%), forests and wooded areas (5.8%), water bodies (3.0%) and built-up areas (1.1%) (based on the BDOT10k topo reference database). As already mentioned, the selection of research areas was preceded by an analysis of the water erosion risk in these areas (Figures 2b,3b,4b).

DTM Source Data
At present, a wide range of DTM and digital surface model (DSM) data is available nationwide (homogeneous data for the territory of Poland). These data were created as a result of the use of various technologies; therefore, they are characterised by varying accuracies and details of mapping the terrain forms, varying spatial resolutions resulting from these parameters, as well as varying validities. The list of the most frequently used data sources is presented in Table 2. As a result of the review and analysis of the above-described sources of elevation data, the following data were selected to test the sensitivity of the soil water erosion model: SRTM, DTED2, LPIS and ISOK.
The DTED2 DTM was developed by the vectorisation of contour lines from military topographic map diapositives with the scale of 1:50,000. The spatial resolution of the raster data was 25 m/30 m. Many tests performed showed that these data are usually of better quality than indicated by the nominal data metric (nominal absolute error in determining the height of points can be up to 26 m). In experimental work carried out at the Warsaw University of Technology, it was estimated that the accuracy of the height determination for the areas around Warsaw is 2 to 3 m.
The LPIS DTM is mainly used to generate orthophotomaps for ARMA (Agency for Restructuring and Modernisation of Agriculture). Aerial photos in the scale of 1:26,000 and topographic maps of 1:10,000 for inaccessible (mainly forest) areas were used for its construction, and new data was obtained in accordance with the data update ARMA schedule. This model has a height accuracy of 1 m-1.5 m. It is possible to obtain free-ofcharge data from The National Geodesy and Cartographic Resources (PZGiK). It is recommended to generate raster DTM with 10-m spatial resolution from the LPIS source data.
In 2010, the implementation of the ISOK project aimed at preventing crises in flood zones was started in Poland. As part of the project, point clouds in the technology of airborne laser scanning (ALS) were obtained. The source measurement data are stored in the form of a classified point clouds with a density depending on the degree of land use from 4 points/m 2 (standard I) to 12 points/m 2 (standard II). As a rule, the data in standard I should meet the accuracy assumptions of the mean error of the x and y coordinates: mx, y <0.5 m and the z coordinate: mz <0.15 m, while the data in standard II should meet the accuracy assumptions of the mean error of the x and y coordinates: mx, y < 0.4 m and the z coordinate: mz <0.10 m [101]. Measurement data, as well as the raster forms of DTM and DSM, can be obtained free-of-charge from the PZGiK resource. It is recommended to generate raster DTM from the ISOK source data with 0.5-m or 1-m spatial resolution, depending on the standard (II or I). This resolution is sufficient to reflect the relief of the terrain with sufficient detail for hydrological studies [90].
The elevation source data was selected so that, in the raster form, they present the sequence of spatial resolution: from 90 m to 1 m (resulting from the accuracy and detail of the source data), as well as different techniques of data acquisition and different ones upto-date. The selected DTMs were acquired on the following dates: SRTM (2000), DTED2 (1980), LPIS (05.2008-06.2009) and ISOK (10.2011-04.2012). Due to the fact that the analysis carried out in this project concerned agriculture areas, the problem of the different validities of the data and model types (DTM/DSM) was considered negligible. All DTM data were preprocessed (e.g., sinks filling or raster snapping) for being as suitable as possible for further analysis.

USLE/RUSLE LS Formulas Implementation
As mentioned in the introduction, few authors have attempted to compare the results of water erosion modelling by testing various data sources for individual criteria. The most frequently used methodologies assume either the study of a new data source or introduction of new approaches to determining the value of one of the factors or the entire model by comparing the obtained new results with the results obtained on the basis of previously used data or technologies [38,42,43,[102][103][104][105] or comparing these with the results obtained using the methodology used so far [17,19,30,36,48,55,106].
The examples of assessing the influence of spatial resolutions of the DTM on the values of DTM derivatives, including the values of the L and S USLE/RUSLE factors [58,88,95,96] found in the literature review, assumed performing tests with the use of a single DTM source with increasingly lower resolutions, obtained by degrading the resolution of the original DTM of high spatial resolution as an adopted methodology. On the other hand, this project assumed the use of various original data sources available for the territory of Poland, which differ in both spatial resolution (resulting from the accuracy and detail of the source data) and the technology of the data acquisition. Initial data processing showed that, despite the similarity of the elevation values, the values of DTM derivatives differ significantly between individual datasets. Therefore, the presented methodology includes not only an analysis of the value of the L and S factors obtained from different source data but, also, a comparative analysis of the DTM values and the values of the land slopes. The first step (stage 1) was to determine the mask of agricultural areas within the boundaries of the research areas. Although the calculation of the values of the L and S factors was performed for the entire studied areas, the comparison of the obtained values was carried out only within the area of agricultural use. This is supported by the fact that water erosion is primarily a threat to arable soils, and the research conducted mainly focused on predicting and preventing erosion in these areas. The second premise results from the previously mentioned specificity of the data used. SRTM data is defined as a digital surface model (DSM), and that is why comparing its derivatives with derivatives obtained from other DTM sources in wooded areas carries the risk of obtaining different processing results. This mask was developed on the basis of data provided by ARMA.
All statistics for DTM data and their derivatives were developed to facilitate a possibility of comparing the obtained results on the basis of the values assigned to points (vector layer) distributed within the study areas in a regular grid of 1m × 1m. Thus, for the results obtained from all DTM sources, the number of observations in each research area is the same. Moreover, this approach guarantees a comparison of the values obtained in exactly the same locations.
The methodology of the comparative analysis of DTM (after preprocessing, e.g., sink filling) and slopes values (stage 2) is presented in Figure 5. The next stage (stage 3) of the study was the calculation of the values of the L factor and the S factor depending on the DTM data source used for the calculations, as well as comparison of the obtained values. An indication of the formulas used and the subsequent analysis included in this stage are presented in Figure 6. When determining the influence of the L and S factors on the water erosion estimation, many authors, referring to the USLE/RUSLE methodology, do not always carry out calculations according to the same assumptions and formulas. The choice of the formulas described by the equations presented in the scheme for the implementation of this stage ( Figure 6) was based on the assumption that the equations should contain as few fixed values and simplifications as possible (e.g., assuming a specific value of one of the parameters for a specific range of slopes). These formulas were also selected by researchers from the Joint Research Centre (JRC)-European Soil Data Centre (ESDAC) [26,29,32,35,41,43,68,92]. The RUSLE2015 project to designate the areas of water erosion risk for EU countries also used them [40,42]. In this study, the simplest single-flow algorithm (D8) was used to calculate the unit slope contribution areas to minimise the impact of the division of the runoff paths or a random factor on the calculation results.
The next stage (stage 4) of the performed analysis consisted of an in-depth analysis of the value of the total impact of the L and S factors, calculated as the multiplication of their values (LS) (Figure 7).
Then, the total LS value threshold for 95% and 99% of the pixels were determined and compared with the maximum values of each of the datasets. The threshold values are understood to be the limit values within which 95% and 99% of the agricultural land in each research area falls, respectively. The significant difference in the value between the threshold value and the maximum value for research indicates the presence of small areas with very high values of the total LS factor values, interpreted as an indicator of the places most exposed to erosion due to their topography.
In the next step of stage 4, the calculated LS values from DTMs with lower resolutions were compared to the LS value calculated from ISOK. The Pearson correlation coefficients were calculated for all pairs of results and the scatter plots for the least and most correlated pairs of the total LS factors values analysed. At the end of this stage, the following statistical values relating to the area of the entire agricultural plot were calculated: the sum, standard deviation, mean, median and the maximum of the total LS factors value. The analysis that complements the described calculations is the visualisation of the DTM elevation value, land slopes and the total LS factors values along the selected terrain profiles. Localisations of the profiles were selected so that they pass through areas with a significant diversification of reliefs and slopes. The comparison of the profiles showed to what extent the DTM data and their derivatives (slopes and total LS factor values) exhibited similarities or differences inside the 3 research areas.

DTM Elevation and Slopes Comparison
A statistical analysis of the elevations in the individual research areas shows that there is a significant similarity of the elevation values in the analysed DTM obtained from four sources. The basic statistics of the DTM values for the research areas are presented in Table 3. Despite the significant differences in accuracy and detail, the comparative analysis of the DTM data shows that the distribution of elevation values is very similar. This is evidenced by the similarity of the histograms and by the high values of the Pearson correlation coefficients. Those values calculated for all DTM pairs are surprisingly high, with (2) the greatest similarity in all research areas shown by the ISOK and LPIS data-the value of the correlation coefficient does not fall below 0.99. The lowest correlation was noted for the DTED2 and SRTM data, but still, the values are not lower than 0.93 (Table 4). Such high values of the Pearson coefficients are also confirmed by the analysis of the scatter plots for the DTM pairs; the closest agreement was obtained for the ISOK and LPIS data, but for all the compared DTM pairs, the points representing the pixel values of the two datasets are arranged closely along the regression line. Despite the very high values of the correlation coefficients, the elevation differences (∆z) of LPIS, DTED2 and SRTM compared to ISOK are significant. For the less-correlated pair (ISOK and SRTM), they reached a maximum of -34.5 m (the elevation of the ISOK is smaller) and 36.3 m (the elevation of the ISOK is higher) for the agricultural areas. Even a comparison of the DTM values for the most correlated ISOK and LPIS pair shows differences on average of ± 0.9 m in all research areas (Table 5). Table 5. Comparison of the differences in the elevation values from DTM of lower accuracy and detail (SRTM, DTED2 and LPIS) to ISOK (minimum, maximum and mean of the absolute value |∆z| of the differences). A statistical analysis of the slope values in individual research shows that there is much less similarity in the value of this DTM derivative than in the case of the elevation values in the four analysed DTM sources. Basic statistics ( Table 6) clearly shows that the maximum slope values increase with the increase of the DTM accuracy and detail (understood as the analysis of changes in the DTM derivative values for the series: SRTM → DTED → LPIS → ISOK). This is a significant change, with the maximum values increasing more than six times between SRTM and ISOK. At the same time, the average value of the slopes and the median grow-this statistical measure is subject to less fluctuations than the average, which means that, although the maximum value increases, high values of landfalls do not occur in a large part of the study area. The greater variability of slope values in the research areas is also evidenced by the increase in the value of the standard deviation along with the increase in the accuracy and detail of the analysed DTM data. The values of the Pearson correlation coefficients of the slope values decrease in comparison to the values of the coefficients calculated for the terrain elevations (DTM values). The slopes calculated on the basis of ISOK and LPIS show the greatest similarity in all research areas-the values of the correlation coefficients do not fall below 0.75. The lowest correlation was recorded for the slopes in the ISOK and SRTM pair; the coefficient values are already lower and range from 0.36 (Bolcie study area) to 0.52 (Janowice study area) ( Table 7). An analysis of the histograms of the terrain slope values also confirms the significant variability in the distribution of the values, with an increase in the accuracy and detail of the DTM. The presence of a few pixels with relatively high slope values in each dataset was noticed ( Table 8). The greater the accuracy and detail of the source data, the higher the maximum slope values. The maximum differences in the slope values (∆s) are significant for all DTM pairs. The largest positive differences (the slope of the terrain calculated on the basis of ISOK data has a greater value) are analogically high for all the compared DTMs and range, on average, from 55° for ISOK and LPIS pair to 60° for ISOK and SRTM pair. Unexpectedly, the greatest negative differences (the slope of the area calculated on the basis of the data from ISOK has a lower value) were obtained for the ISOK and LPIS pair: -40.4° for the Janowice research area. The mean values of the absolute differences (|∆s|) do not differ significantly between the individual DTM pairs, ranging from ±2.6° for the ISOK and LPIS pair to ±4.6° for the ISOK and SRTM pair (Table 9).

ISOK-SRTM (m) ISOK-DTED2 (m) ISOK-LPIS (m)
When analysing the distribution of the ∆s values, the greatest values are observed in the middle and lower parts of steepest slopes. Panagos [40] mentioned that the vast majority of areas in Europe have slope values of up to 50% (approx. 26.6°). This is confirmed by the results obtained for the research areas, but only if the results of the analysis using DTM data of similar accuracy and detail are compared (the analyses described in the article were carried out using EU-DEM with a spatial resolution of 25 m). The values obtained for the research areas for the DTED2 slope maps reach 29°, while, for DTM with higher accuracy and detail, the maximum values for the Janowice area (foothills) reach 46° (for LPIS DTM) and 78° (for ISOK DTM). Table 9. Comparison of differences in the values of the terrain slopes (∆s) (°) calculated on the basis of DTM (SRTM, DTED2 and LPIS) against ISOK (minimum, maximum and the mean of the absolute value |∆s| of the differences) in the 3 research areas.

Comparison of the Total LS Factors Values
The specificity of the universal soil loss equation USLE/RUSLE results in the fact that the influence of the variability of the values of the L and S factors on the value of soil loss due to water erosion is multiplied. The total influence of the L and S factors in the soil loss equation is determined as the multiplication of their values (L * S). The first part of comparing the total LS factors values was to compare the results for each pixel (1m x 1m grid). The second part was to compare the statistical values after aggregating for the agricultural plots. For individual research areas, the maximum LS range was obtained up to the value of 3314.9 (for the Janowice research area) to 824.3 (for the Winiarki research area) and to 1116.8 (for the Bolcie research area); the higher values were obtained based on ISOK DTM data (Table 10). In all the research areas, the maximum values of the total LS factor values increase consistently with the increasing accuracy and detail of the DTM source data. In the research area of Janowice, the mean and median values, as well as standard deviation, clearly decrease, along with the increase in the accuracy and detail of the source data.
There is no such relationship in the Winiarki and Bolcie research areas, and the mean values and standard deviation are the highest for the LS calculated from the DTED2 data. This is confirmed by the analysis of the terrain profiles presented in Figure 8 and in Appendix B. A comparative analysis of the elevation profiles indicates that, along with a decrease in the accuracy and detail of DTM, a process of flattening of the valleys takes place, which leads to a decrease in the length and steepness of the slope represented by the generalised source data compared to the high-resolution data. Along with the increase in the accuracy and detail of the DTM, a significant increase in the area with the total LS factor values close to 0 is noticed. This is evidenced by the total LS factor threshold values obtained for 95% and 99% of the area. The highest threshold values in the Janowice research area were observed for SRTM derivatives and, in the Winiarki and Bolcie research areas, for the DTED2 derivatives (Table 11). There is a significant discrepancy between the total LS factors threshold values for 95% and 99% of the area and the maximum of the total LS values, and this discrepancy increases with the increasing accuracy and detail of the source DTM.
The analysis of the spatial distribution and the histograms of the total LS factor values shows that the largest range of the areas with the lowest values was obtained for the results calculated on the basis of the ISOK data ( Figure 9). Areas with values ranging from 0 to 1 occupy an average of 60% of the research areas for the results obtained from SRTM and 88% of the research areas for the results obtained from the ISOK data.  The analysis of the values of the Pearson correlation coefficients (Table 12) and the scatter plots indicates a significant lack of similarity of the total LS factors values obtained for the pairs of different DTM sources. Still, the lowest correlation values were obtained for the ISOK and SRTM pair-the Pearson coefficient value does not exceed 0.08 in all the research areas and indicates a complete lack of similarity of these datasets. The maximum value of the correlation coefficient was obtained for the DTED2 and SRTM pair for the Bolcie and Winiarki research areas, but it is still only a maximum of 0.44.  The distribution of the values in the profiles showing the total LS factors confirms the previously described discrepancy in the values calculated on the basis of different DTMs. In most cases, it is difficult to find the convergence of the local maximums and minimums in the graphs. There is also a noticeable, though less visible in the previously analysed statistics, overestimation of the total LS values obtained on the basis of the LPIS data-much more noticeable in the research areas of Winiarki and Bolcie. The profiles of all the compared features (elevation, slopes and the total of the LS factors values) for each of the research areas are presented in Figure 8.

Comparison of the Total Values of the LS Factors in Agricultural Plots
The next stage was an analysis of the total LS factor values after aggregating the results in the agricultural areas (plots) designated on the basis of the data provided by ARMA. For all the agricultural plots, the LS statistical values were calculated relating to the area of the whole plot: sum, standard deviation, mean, median and maximum values. These values are presented in Table 13. Visualisation of the changes in the value of the sum of the total LS factor values in the plot and the maximum value of the total LS factor values in the plot are presented in Figures 10 and 11. Visualisations of the total LS factor values for the individual plots ( Figure 10 and Figure 11) in all the research areas show: • a decrease in the value of the sum of the total LS factor values in the plot, • a decrease in the percentage of plots with high values of the sum of the total LS factor values and • a significant increase in the maximum of the total LS factor values in the plot as the accuracy and detail of the source DTM increases.
The obtained values of the Pearson correlation coefficients for the plots indicate a significant increase in the similarity of the total LS factor values analysed in the areas of the agricultural plots in relation to the analysis of the values at the points placed in the regular grid (1m x 1m). The highest agreement was obtained for the sum of the total LS factors values in the plot for the ISOK and LPIS pair; in all the research areas, it is a very high similarity-the values of the correlation coefficients have a value of approx. 0.96. Even for the pair of ISOK and SRTM, although the lowest, they range from 0.52 to 0.78 (Table 14).   Significantly lower values of the Pearson correlation coefficients were obtained when comparing the maximum values in the agricultural plots. These values significantly differ from each other, and the highest correlation occurs for the DTED2 and SRTM pair at the level of 0.43-0.63 and the lowest for the ISOK and SRTM pair at the level of 0.07-0.49 ( Table 15). The difference in the total LS factors values in an agricultural plot is indicated by the results of the variation index Wz, which was used to compare the LS values obtained from the SRTM, DTED2 and LPIS data with the LS values obtained from the ISOK data (Equation(2)).
As already mentioned in the description in the Methodology section, very low values of the Wz index (close to 0) indicate a significant underestimation of the total LS factor values relative to ISOK; in turn, very high Wz values indicate a significant overestimation of the total LS factors values relative to ISOK. The values of the similarity index ranging from 0.9 to 1.1 (±10%) were considered as values comparable to those obtained from ISOK. The comparison of the sums of the total LS factor values in the agriculture plots calculated on the basis of SRTM, DTED2 and LPIS against the sums of the total LS factors values calculated on the basis of ISOK for the three research areas is presented in Table 16. The calculated values of the coefficient of variation Wz range from 0.001 (the total LS factor values with SRTM were 1000 times lower than the value obtained from ISOK) to 475 (the total LS factor values from SRTM were 475 times higher than the value obtained from ISOK). Both extreme values were obtained for the Janowice research area. The greatest differences were noticed when comparing the total LS factor values obtained from SRTM with the values obtained from ISOK. The analysis of these values revealed a significant percentage of the areas with an underestimation the total LS factors values; in all the research areas, it is a significant part of the area (Table 16). This relationship is also visible in the visualisations of the distribution of the WZ variability index ( Figure 12) calculated for the total LS factors values with DTED2 and ISOK, as well as LPIS and ISOK.

Discussion
Despite the advantages of the models based on the USLE/RUSLE methodology, one should also be aware of their limitations. Many authors point to the mismatch between the model and the European conditions. The empirical model was developed for the eastern part of North America, which, compared to Europe, is characterised by a different distribution and nature of precipitation, a different landscape and a different distribution of the dominant hydrological processes. Moreover, the model does not take into account gully erosion (ephemeral gullies) and crop erosion [10,40,107]. Some researchers report that comparing the modelling results with the field studies indicates that the USLE/RUSLE models usually give slightly overestimated values of soil losses or there is a noticeable overestimation of low values and an underestimation of high values [32,43,69,84]. Another important problem and challenge in erosion modelling is the assessment of the accuracy of the obtained results. It is a very difficult and complicated task, as it requires an estimation of the uncertainty of each of the equation factors, including the L and S factor values [86]. It was also noticed that the modelling results are sometimes inadequate for long slopes of complicated shapes [107]. The use of a more detailed DTM increases the length of the slopes and the complexity of their shapes. For a more complete analysis of the impact of DTM sources on the modelling results, the results obtained with the use of other, multiple-flow runoff algorithms should also be compared, as well as eroded soil deposition areas, as indicated by the works [37,60].
The research carried out indicates that the overestimation of total LS value calculated for the DTM with lower accuracy and detail concerns a much larger area of land than the underestimation that occurs in relatively small areas. However, the values of the differences obtained in this case are much higher. This is probably related to the less accurate determination of the water paths flowing down the slopes and unit upslope contributing areas using a DTM with less accuracy and detail. The obtained conclusions confirm the observations of Evans [107], who compared the results of field studies with the results of modelling using a DTM with low accuracy and detail. He showed that modelling using the USLE/RUSLE methodology results in the overestimation of values where low total LS factor values were found in the field and an underestimation of values where high total LS factor values were found in the field. Moreover, he notes that significantly high soil losses occur in small fragments of farmland. These observations are confirmed by the modelling results obtained in this study. With too-low accuracy and detail of the DTM, high values of the L and S factors are obtained, while, with a DTM with high accuracy and detail, low values of L and S are actually obtained in most of the area and very high values in small areas of the field, which indicated the areas of accumulation of runoff on slopes with significant steepness. Similar conclusions were drawn by Yang et al. [108], Fu et al. [109] and Raj et al. [110]; however, the above-mentioned studies were conducted outside Europe.
The comparison of the values of the LS factors indicates that the values obtained from the source DTMs with different accuracies and details and, resulting from them, different spatial resolutions, do not change uniformly with the increase in the DTM accuracy and detail. These dependencies cannot be described in a simple way, because for the same region, both high total LS factor values from one source DTM are obtained for the low total LS factor values from another DTM source and vice versa.
Opinions on the overestimation of the values of soil loss in the USLE/RUSLE models resulting from the literature review may be related to the inability to include natural barriers in the terrain in the available algorithms (roads, paths, fences, bushes, etc.), which are natural blockades for runoff water, often causing eroded material to be deposited. This was also mentioned by Panagos [42]. However, it seems that DTM data from airborne laser scanning (ALS) partially solves this problem.
Thus, there is a need to calibrate the existing erosion empirical models under the new conditions of access to data with high and very high accuracy and detail, as mentioned by Mitosova et al. [84].
In the literature on soil water erosion modelling, opinions can be found that the values obtained from modelling should be treated as a potential risk (qualitative approach) rather than the actual values of soil losses (quantitative approach) [12,52].
The analysis of the distribution of the obtained results of the differences in the values did not allow for a certain formulation of conclusions about which locations these differences are the greatest (e.g., the value of the slope and location in a specific part of the slope).
For the presented study, the ALS data was obtained with a scanning density of four points/m 2 , so generating a DTM with a resolution of 1m x 1m resulted in additional averaging of the elevation values and partially reduced the possible noise. Further data processing was not carried out on purpose, so as not to change the source data values. The next direction of research could concern ALS data filtering and comparing the results in the field studies, both the DTM values and the slope values, as well to correctly define the surface runoff paths.

Conclusions
As a result of the analyses carried out in accordance with the presented methodology, a significant impact of the accuracy and detail of the source elevation data on the results of calculating the values of the factors of the universal soil losses equation-slope length (L) and slope steepness (S), considered as the key of the USLE/RUSLE model factors-was demonstrated.
A different image of the terrain surface presented with the use of various DTM models (SRTM, DTED2, LPIS and ISOK), despite the high correlation of the terrain elevation values, shows a significant differentiation of the DTM derivative values. It has been proven that, when using the DTM with the lowest accuracy and detail (e.g., SRTM)-often recommended for analyses for the entire country (in accordance with the principle of "from general to detail")-not only a significant (even a thousandfold) underestimation in one part of the studied area could be expected but, also, an overestimation in another part of the studied area of the LS factor values.
The conducted research has shown that the study of local phenomena (and soil erosion is one of them), even when carrying out an analysis over a large area, should use data of the highest accuracy and detail, because the use of data with a lower level of detail and accuracy may lead to wrong decisions. Therefore, it seems to be wrong to think that the detail of the source spatial data for an analysis should be selected adequately to the scale of the target study area, i.e., for example, in studies of the country area, data should be used at a certain level of generalisation. The modelling of local phenomena should use data adjusted to the scale of the phenomenon (and not to the scale of the study), and the visualisation of the results in smaller scales can be developed using appropriate cartographic generalisation methods.
The ALS data are characterised by high scanning density (in Poland, data from the national resource have a standard density of four or 12 points/m 2 ). More research is needed to determine the optimal scanning density, data filtration and final raster data resolution for the proper modelling of surface water runoff phenomena.
It is necessary to pursue a responsible land use policy using the advantages of GIS technology. The identification of vulnerable areas should be performed using data of adequate accuracy and detail to prevent negative changes and protect those sites where there is a real risk of soil water erosion, and it is undoubtedly necessary. It is not possible with the use of source data that do not reflect the topography as accurately as possible. It is also important not to limit agricultural activities and not to impose special-other than those commonly used in a given area-treatments, if it is not necessary. These decisions also have financial consequences, both for agricultural producers and for the state budget.
Despite many years of research and technology development, there are still many difficulties in recognising, describing and quantifying soil water erosion, and there are still limited possibilities for predicting the intensity and frequency of erosion events [12,43,87]. In a comparative analysis of the obtained results, it may be helpful to develop a set of metadata, the provision of which would clearly indicate the source data characterisation, algorithms used and the way of their implementation. The proposed set of such metadata is included in Table 1. Moreover, another difficulty points to the up-to-date conducted research. Rainfall, land cover and topography data acquisition periods are often different, and some anti-erosion treatments may change from year to year.
Thus, some questions Boardman asked in 2006 [87]: "Where is erosion occurring? Why is it happening? How serious is it? Who does it affect? What should be the response? Can we prevent it?" are still relevant.
Funding: This research received no external funding.

Acknowledgments:
The author acknowledges the Agency for Restructuring and Modernisation of Agriculture (ARiMR) for providing the plot data and the Institute of Soil Science and Plant Cultivation (IUNG) for providing the expert erosion risk estimation data used in the analysis. I would like to take this opportunity to acknowledge the time and effort devoted by the journal reviewers to improving the quality of my published work. I appreciate the detailed and useful comments that I obtained from them.

Conflicts of Interest:
The author declares no conflicts of interest.

Appendix B
Visualisation of the elevation (m), slopes (°) and the total LS factor values (Y-axis) calculated from all the DTM sources for a selected section of the profile (over a 400-m section) for the Janowice and Bolcie research areas (visualisation for the Winiarki research area is presented in Figure 8).