Development of Topographic Factor Modeling for Application in Soil Erosion Models

1. Introduction

The fast changes in soil use and higher vegetal resource demands favor the triggering of water erosion that needs to have its rates expressed in space and time for the proper adaptation of control practices and resources for agricultural planning. The physical processes of disaggregation, transport and soil deposition that define the erosive process are hydrologically directed and the movement of the water on the soil undergoes the interference of the topography, climate, soil class and land use, so that the studies regarding the theme are based on the intense experimentation of the effects of the variations of these factors on the sediment production. The estimate of the topographic variables, although benefitted by automatic generation and spatial distribution made possible by the Geographical Information Systems (GIS’s), is the target of controversy related to the formulation of algorithms for this end, so that its three-dimensional calculation is not a current procedure in the geoprocessing programs [1].

Given the scope of the process, the topographic modeling in the erosion analyses can differ in terms of complexity, processes considered and data required for model use and calibration, which can be empirical, physical and conceptual [2]. In the empirical models most used, such as USLE (Universal Soil Loss Equation) [3] and the revised version of USLE (Revised Universal Soil Loss Equation - RUSLE) [4], the topographic factor is expressed by the association of the steepness and the length of the slope called, respectively, factors S and L. Considering the proper formulation of USLE and its adaptation to the work context in the Digital Elevation Model (DEM), obviously the advantages associated to DEM derive almost entirely on issues related to the topographic LS factor, because it can be evaluated with the aid of DEM and where the precision of the extracted parameters can become apparent [5].

An aspect that hinders the estimate of appropriate topographic factor (LS) values for applications in GIS and results in high limitation in the use of the USLE and RUSLE erosion models [6, 7] are the dynamics of the erosive process in complex reliefs and hydrographic basins, since USLE was primarily developed for the prediction of the erosion in not very accentuated and uniform slope stretches, in other words, not considering if they are concave, convex [8], or in combination. The limitation of the empirical modeling in the soil loss estimates in complex profiles impelled the development of conceptual models (or semi-empirical), such RUSLE 3D and USPED (Unit Stream Power Erosion and Deposition) [9]. Derived from USLE, these models intend to represent their advancements by adding a physical basis that tries to relate the morphology of the relief and the erosion defining parameters.

In this sense, given the strategic need for generation and diffusion of algorithms for automatic mapping of the topographic variables used in the operationalization of the digital analyses of water erosion, the objective of the present chapter was to conduct a review of the topographic factor development in erosion equations applied in computational geoprocessing systems, with prominence for USLE and RUSLE, addressing the main theories and algorithms used in the digital treatment of the data.

2. Digital Elevation Model (DEM)

In the landscape, the topography determines the behavior of the surface runoff, the phase of the hydrologic cycle that is most directly associated to the water erosion and that requires a rigorous and effective analysis throughout its entire extension, make possible with the use of digital elevation models (DEM). The analyses developed on a DEM allow: to visualize the model in planar geometric projection; generating gray scale images, shaded images and thematic images; calculating fill (embankment) and cut volumes; conducting profile analyses on predetermined trajectories; and generating derivative maps, such as steepness and exposure maps, drainage maps, contour maps and visibility maps. Products of the analyses can even be integrated with other geographical data types aiming at the development of several geoprocessing applications, such as urban and rural planning, agricultural suitability analyses, risk area determination, environmental impact report generation [10], elaboration of digital soil maps, as well as maps of soil attributes such as soil organic matter content [11], among others. Therefore, DEM should faithfully represent the relief allowing to capture the topographic variations presented.

The elaboration and creation of a DEM, indispensable for the representation of a real surface on the computer, can be represented by analytical equations or a network (grid) of points, in a way that transmits the spatial characteristics of the land to the user [12]. Therefore, the information contained before in specific points (vectors) are transformed into a continuous spatial distribution of the relief (raster), enabling new inferences about the local relief. Different methods exist for the interpolation of the data and DEM generation, which are built through regular rectangular grids, such as the Topogrid [13], or triangulated irregular networks (TIN) [14]. For the choice of efficient DEM in the evaluation of the erosive process, an intense preliminary analysis of information, from a hydrologic point of view, is recommended, because the development of the water erosion occurs in response to the manner the water moves through and on the landscape [15].

The geomorphological and hydrological consistency of a DEM is reached when the matrix image faithfully represents the relief features, such as the hydrographic basin watershed, thalwegs and concave and convex elements, and it assures the convergence of the surface runoff for the mapped drainage network. In this sense, several water erosion analysis and modeling works have used the TIN model [16-18], as well as the Topogrid model [19-21] for DEM generation. In a research conducted with the objective of defining the drainage network in a sub-basin [22], the original contour curves (scale 1:10.000) were compared to the curves generated by DEM's of the Topogrid, linear TIN and TIN natural neighbor interpolators. A higher Topogrid hydrological consistency was observed verified in the better continuity of the contour curves and higher drainage area and watershed detailing, resulting in a smaller amount of flat areas and in more detailed drainage pathways. The authors emphasized that all of the generated models have high reliability due to the precise topographic surveys data from which they originated, a fact also observed by [15, 18]. Starting from this same comparison among interpolators, other works [23, 24] made similar observations regarding the behavior and reliability of the models.

The precision of the data collection will influence the quality of the corresponding digital model and the choice of the database to be used in the construction of DEM becomes fundamental. Such data can be derived from contour curves, elevation points, photogrammetric analysis from aerial photography, information collected by stereoscopic satellite images, or radar [25, 26]. The digital database can result from a digitized manual survey obtained through direct readout digital equipment (total stations, topographic GPS), or by remote sensing equipment (radar, laser) in which different DEM generation models can be applied supplying DEM's with varying precision [26]. The spatial distribution and the amount of errors propagated can vary according to the spatial resolution, therefore, the effect of the spatial resolution on supplying useful information for the determination of an appropriate resolution should be investigated [27]. The higher the scale of a map, the higher detail it presents, and to the contrary, the lower the scale in a map, the higher the degree of generalization seen which increases the minimum cartographic area of the land [27]. With better horizontal resolution DEM’s, it is possible to include more relief roughness aspects, reducing the length of the slope straight line segments and increasing the accuracy of the L and S factors [18].

The evolution of the geographical information systems and the growing availability of better quality radar images enabled the obtaining of terrain elevation models (DEM) with increasingly better spatial resolutions. However, data surveys considered as having high precision (below 10 meters of resolution) still possess high costs, limiting their use in research. Currently, the most widely used digital database for the generation of DEM’s originates from of the digitization of topographic maps or those obtained by remote sensing maintaining their original precision. Because the influence of the pixel size has a significant weight in the analyses derived from DEM’s, the choice of the spatial resolution proportional to the scale of the primary data must have certain considerations, among them, the original contour curve scale and the characteristics of the mapped relief. For instance, with a minimum horizontal distance between the curves on the order of 20 m, the spatial resolution of 15 m can be shown appropriate for the detailing of the relief presented in the original base, considering that lower resolutions would tend to generate erroneous information (nonexistent), while higher resolutions would not detail the relief in a satisfactory way [11]. Furthermore, a sufficient spatial resolution cannot only depend on the aimed information and/or the precision used in the collection of this information, but also on the topography. In areas where simple hillsides exist with flat topography, a coarse resolution (> 20 m) may not lead to major errors in hydrographic basins, being able to be used with little uncertainty. In a complex topography, with accentuated slopes, a coarse resolution can result in great uncertainty, and a better resolution could be necessary.

Along those lines, studies have been developed seeking to define the best DEM resolution that is capable of precisely representing the variations of the relief, thus reducing the uncertainties of the erosion prediction models that need topographic modeling. In Slovakia, the S factor derived from a MED with 50 m resolution, obtained from digitized contours of the topographic maps in a 1:50.000 scale, promoted a sufficient level of detail for this type of regional evaluation and the spatial resolution selected reflected the scale of the primary data [28]. To determine the resolution of DEM when the objective is to promote an impartial average global estimate, the global variance of the LS factor can be used [29].

The local variance of a cell measure the average local space variability, while the global variance measures the global variation of the estimates showing a considerable difference of behavior with the increase of the DEM resolution. Thus, in [29] the objective of the research was to evaluate the appropriate DEM resolution for the spatial prediction of the LS factor. A decrease of the global variation of the LS values obtained in the 30 by 100 meter resolution in DEM was observed, followed by a stabilization after 100 m. The local average variance and semivariance of a cell increased with the 30 by 50 m resolution and it decreased after 50 m. The highest local variance and semivariance of a cell was at the 50 m resolution and thus it could be considered appropriate for the necessary detailing of the spatial information (distribution and variability) for the LS factor prediction [29].

A study developed in Thailand analyzed the influence of the spatial resolution on the results of the LS factor [30]. Two DEM resolutions extracted from a SRTM (Shuttle Radar Topography Mission) radar image with 90 m of resolution (original image resolution) and 30 m (resampling of the 90 m resolution) were appraised. The grid size change affected the steepness values, compromising the L and S factor values, since the L factor depends on the grid size and the steepness and the S factor only on the steepness. When affecting the L and S factors, the resolution also affected the sediment transport ratio. The best sediment production estimates were observed in DEM with resolution of 30 m. A fundamental observation is made by the authors of the study, who highlight that the better results of the 30 m resolution compared to the 90 m using the USLE methodology, is probably due to fact this resolution is closer to the 22.4 m slope length, the length used in the derivation of the USLE relationships. Another study [31] compared SRTM radar images with spatial a resolution of 90 x 90 m and digitized hypsometric curves to determine the topographic factor. The highest detail was obtained by the SRTM images that evidenced lower LS ranges. They concluded that the difference occurred due to the higher detailing in plane areas, where LS is lower (ramp height lower than 40 m, with low steepness).

3. Water erosion modeling

3.1. USLE and RUSLE empirical modeling

The empirical models are the simplest ones, generally possessing less data and lower computational base than the physical and conceptual ones. The empirical models normally have a high aggregation of time and space and are based on analyses of the erosion process using statistic techniques. For this reason, they are particularly useful as the first step to identify the sediment sources. Universal Soil Loss Equation (USLE) and Revised Universal Soil Loss Equation (RUSLE) have been the most used models in the world for predicting erosion processes due to their simplicity and the availability of information.

The topographic factor is the most sensitive parameter of USLE/RUSLE in the soil loss predictions, where a higher relative effect of the steepness factor is observed in a simple analysis of sensitivity. However, an interaction of the steepness (S) and the slope length (L) exists and, LS being used as the “only” parameter in the sensitivity analysis, its influence is even higher on the soil loss than the remaining parameters, including L and S individually [32]. By definition, the slope length (L) is the distance from the point of origin of the surface flow to the point where each slope gradient (S) decreases enough for the beginning of deposition or when the flow comes to concentrate in a defined channel [3]. The soil losses increase with the increase of the slope length and steepness, conditions where the surface flow reaches high-speeds.

The procedure to obtain the slope length was originally manual in these models, which may be adapted to GIS framework. It initially consists of slope length identification through information plans of slope steepness and classified aspects, such as the exposition angle between the slopes and the north. From each rill slope or polygon, the average slope steepness (degrees) and the altitude are calculated [33, 34]. It is possible to calculate the slope length through the following equation [33]:

$$\text{L}=\frac{\text{DH}}{\sin \alpha }$$E1

Where L is slope length (m); DH is altitude difference (m); and α is average slope steepness (degrees).

The slope angle α corresponds to the inverse of the tangent angle that may be calculated dividing sin α (altitude difference) by cos α (distance between level curves and/or quoted points). The α angle of a surface defined by two points (A and B) is calculated with the horizontal as show in Figure 1 [35].

Figure 1.

Trigonometric variables in the calculation of the slope.

The slope length (L) and slope steepness (S) factors of the USLE were developed for uniform slopes based on empirical models, which means that they use dependent field measurements. For USLE/RUSLE, they are calculated from the comparison with a ramp length of 22.1 m and 9% slope with the use of a factor m for different steepness classes [3].

The L factor can also be obtained by pixel size in the DEM. If one size of the pixel is considerate as flow length, an equal value is determined for the flown length; therefore it is assumed that the inclination is composed by segments of equal dimensions, consequently with different inclinations, which is not true. However, this approach is considerably feasible if a pixel of suitable dimension is used [36]. A study developed in Thailand evaluated two DEM resolutions and observed better results from 30 m resolution using USLE methodology due to the fact that this one is closer to 22.1 m slope length, which is used for the derivation of model relations [30].

The calculation methodology of LS factors proposed by the USLE was improved in the equation revision, named RUSLE [4], considered more extensive than the previous model. The L factor, in both USLE and RUSLE, is expressed as [3]:

$$\text{L}={\left(\frac{\lambda }{22.1}\right)}^m$$E2

Where L is slope length effect on soil loss standardized for 22.1 m length; λ is field slope length (m); and mis slope length exponent.

In the USLE, the m recommend value is from 0.2 to 0.5 for slope levels lower than 1%; 1-3%; 3.5-4.5%; and 5% or more, respectively. Therefore, if a slope gradient is higher than 5%, slope length factor do not change with slope inclination. However, in the RUSLE, m continues to increase with slope inclination (Equation 3). Thus, in the RUSLE, the slope length effect is a function of the erosion ratio of rill to interrill [36].

$$m=\frac{\beta }{(1+\beta )}$$E3

$$\beta =\frac{\left(\frac{\sin \theta }{0.0896}\right)}{\left[3{(\sin \theta )}^{0.8}+0.56\right]}$$E4

Where β is ratio of rill to interrill erosion; and θ is slope angle.

Researchers observed that the m = 0.5 exponent of USLE is better adapted for very accentuate slopes [37]. When the slope increases from 9% to 60%, the mexponent increases from 0.5 to 0.71. The slope length exponent, m, is 0.7 for a 50% slope with 60 m length and a more moderated ratio of rill and interrill erosion. When the 0.7 factor is used, the RUSLE predicts an addition of 22% of soil loss than the USLE (m = 0.5) through a 60 m of length slope. When the slope is lower than 9%, the USLE will predict a higher soil loss than RUSLE and, being steeper than 9%, the RUSLE will predict a higher soil loss than USLE. The higher difference occurs in much accentuated slopes.

The equation used in the USLE slope factor (S) (Equation 5) was modified to obtain more accurate results in RUSLE model (Equation 6), probably due to changes in the slope factor [34, 38], which depends on the slope angle Ө.

$$S=65.4{\sin }^2\theta +4.56\sin \theta +0.0654$$E5

$$S=10.8\sin \theta +0.03,\text{for}\theta \text{<}\text{9\%}\text{or}S=16.8\sin \theta -0.50$$E6

3.1.1. USLE and RUSLE limitations

The restrictions of the USLE and RUSLE empirical models frequently occur because neither examines the hydrologic phenomena in their geographical context, using a simplified representation of spatial elements that assumes the hydrographic basin as uniform [39]. Many methods have been developed seeking to include complex slopes, common in a context of hydrographic basins [40]. In a comparison of several manual methods it was concluded that there is no obviously better method [41]. The errors of the empirical models are produced because the water erosion, being a hydrologically driven process, is not evaluated in relation to the surface runoff [42-45]. In the soil loss estimates using the USLE and RUSLE models the surface runoff is not considered in a direct way, though they indirectly consider that the flow transports the eroded sediment and the concentration of sediments depends on the kinetic energy level of the rain, in the sample space of a parcel [44]. Thus, the surface runoff in the empirical models is a primitive factor. This presupposition limits the potential of these models in predicting erosive factor changes, on the scale of basins or drainage systems, which are favored in models based on physical and semi-empirical processes where the surface runoff constitutes a fundamental factor in the water erosion prediction.

For local conservation planning, the LS factor is usually estimated or calculated from length and inclination measurements in the field [6], or even through manual procedures on cartographic bases, making the procedure very difficult and slow due to the difficulty of individualization of each slope [1]. The measurement of the ramp length is made from the evaluated point in relation to the watershed. Besides possessing the difficulty of locating the watershed, this procedure considers the straight line distance until the watershed, concealing the importance of the relief form, because the erosion is affected by the torrent that comes from the whole contribution area. These labor intensive in-field measurements rendered the soil erosion modeling obviously unviable on a regional scale [6] leading to the determination of the ramp length based on the estimate of an average value for hydrographic basins, which is an oversimplification of the true situation [7]. Furthermore, an underestimate of LS values, obtained manually, and consequently also of the erosion risk is observed when compared to the irregular slopes considered in automated models [46].

A second shortcoming of these models is the evaluation of only the erosion, without sediment deposition prediction [47]. When adopting an average rate for an entire slope or hydrographic basin, addressing the erosion using the USLE and RUSLE models does not offer any information as to the sources and sinks of the erosion materials. In spite of the methodology of dividing complex landscapes into series of semi-homogeneous planes used by these models, to provide some consideration as to the convexity and concavity of the inclination, the erosion is only calculated along the flow in a rectilinear manner, without full consideration of the convergence and divergence flow influence [45]. No approach adequately supplies spatially distributed information on the erosion necessary for effective control of the erosion and sediments. Thus, on a hydrographic basin or landscape scale, the spatial distribution of the soil erosion predicted by such models will distort the current conditions and will tend to overestimate the erosion [47, 48]. Some studies mention overestimates in lower soil losses and underestimates for the high losses using the USLE and RUSLE models [44, 49, 50]. As a solution, studies recommend to first identify those portions of the landscape subject to the deposition and to exclude them from the analysis when applying the USLE and RUSLE models [9].

USLE was related to GIS due to the advantages of handling great amounts of spatial data. Until the middle of the 1990’s, a great limitation in the use of the USLE and RUSLE erosion models on a regional landscape scale was the difficulty in estimating appropriate LS factor values for applications in GIS [7], since such models evaluate the effects of the topography on the erosion in a two-dimensional way. In that context, the use of models distributed in space came to represent a powerful environmental analysis tool, highlighting soil erosion by water on the hydrographic basin scale.

3.2. Conceptual modeling

The conceptual methods incorporate the impact of different erosive processes through empirical parameters [51] usually obtained through calibration with observed data, such as flow discharge and sediment concentration [52]. Therefore, these models represent the processes within the scale in which they were simulated [53]. It is noteworthy, particularly on a large scale, to mention that deposition patterns and sediment residence time are still little understood in a way that the erosion prediction and the sediment deposition rates on these scales are based, usually, on empirical or semi-empirical studies that are applied in a uniform way throughout the whole area [54].

The semi-empirical LS factor explains the double phenomenon of drainage convergence and furrow [27]. The result of the LS factor thus comes to be equivalent to the traditional LS factor on flat surfaces, but with the advantage of being applicable to slopes with complex geometries [55-57]. When substituting the empirical topographic factor by the semi-empirical one in USLE, the laminar and concentrated flow in complex terrains is considered in the spatial distribution of the erosion, making the estimate more precise.

In the conceptual models the slope length factor is substituted by the upstream contribution area [9, 46, 55, 56] whose modeling conducted in the digital elevation model (DEM) allows to determine the drainage network considering the direction of the surface runoff and the accumulated flow. For each cell, the contribution area upstream is obtained from DEM initially calculating the steepness and aspect maps, building the water flow paths later. The upstream contribution area map is determined from the water flow path lines and the DEM spatial resolution. The precision of the model is related to the uncertainty of the empirical parameters used in the LS factor equation, the accuracy and resolution of DEM and to the methods for derivation of the topographical variables related to LS, such as steepness, aspect and contribution area [27]. In that way, the topographical LS factor can be finally obtained.

Incorporating this concept, an equation modified to compute the LS factor in the form of finite difference in a grid of cells representing a segment of the hillside was derived [46]. Another model, called RUSLE 3D (Revised Universal Soil Loss Equation 3D) presented a simple and continuous form of the LS factor equation considering the impact of the convergent flow [9]. Also considering the contribution area, the USPED model (Unit Stream Power Erosion and Deposition) was developed from the drainage force unit theory [56, 57] for analysis of the erosion and deposition.

3.2.1. Contribution area modeling

The modeling of the contribution area is conducted resorting to DEM, because it contains information that allows to determine the surface runoff network. As such, based on DEM, the flow direction and the accumulated flow and the steepness are determined. The area of contribution of each cell (pixel) of DEM, considering a grid of cells, is its own area plus the area of the upstream neighbors that possess some drained fraction for the pixel in question. The contribution area (A) of a specific grid of cells is calculated from the product of the accumulated flow (χ) and the area of each cell (η) [58]:

$$A=\chi \eta$$E7

The determination of the accumulated drainage areas (or accumulated flow), which allows the simulation of the hydrographic network, are defined based exclusively on the flow directions. The accumulated flow represents the amount of rain that will drain through each cell, supposing that all of the rain become torrents and there is no interception, evapotranspiration, or loss of underground water. Each pixel receives a value corresponding to the sum of the areas of all of the pixels whose drainage contributed to the analyzed pixel [59].

The flow direction defines the flow path of water, as well as sediments and nutrients, in areas adjacent to the lower altitude points in all of the positions in the hydrographic basin [60]. Independent of the magnitude of the rain event, the flow algorithm in a GIS establishes a one-dimensional flow network connecting each cell with other cells of the hydrographic basin in DEM until the point where the whole surface runoff generated inside the hydrographic basin meets, defined as the mouth [61]. As such, the hydrological relationships are built between different points within a hydrographic basin, topographical continuity being necessary so that functional drainage exists [62].

The estimate of the flow direction is based on the physical principle that the mass of controlled gravity proceed in the direction of the most accentuated slope. The slope is characterized identifying the plane tangent to the topographical surface in the center of the cell. The maximum plane elevation change rate characterizes the inclination gradient, while the correspondent cardinal direction of this larger difference is the aspect [63] (Figure 2).

Figure 2.

Code of flow direction in analogy to cardinal points generated on the aspect map.

Many algorithms have been developed to consider the contribution area. The flow direction methods, profoundly different, are classified in: concentrative, also called single direction or eight directions, that transfer the whole source pixel matter to the downstream pixel in question [59, 64-68]; and dispersive or multiple direction, that divide the matter of the source cell among several receptors [63, 66, 69-73]. That means that, while those of single flow consider all of the slopes as concave or parallel, the multiples also differentiate the convex ones [46]. In water erosion analysis, the most thoroughly used methods are the concentrative [1, 7, 29, 47, 74-81].

The first and simpler method to specify the flow direction attributes the flow of each pixel to one of their eight neighbors, be it adjacent or diagonal, in the direction of the steepest hillside slope. This method, designated Deterministic 8 Algorithm (D8), is based on the fact that the water can move in 8 possible directions, as demonstrated in Figure 2 (8 flow directions) [64]. The D8 approach has disadvantages arising from the determination of the flow distributed equally in only one of the eight possible directions, separate by 45˚, that is expressed in parallel (or convergent) flow patterns in the directions of the cardinal or diagonal points, intermediate values not being possible [60]. In a complex topography, however, the divergent flow frequently can occur causing a significant impact on the delimitation of the basin contribution area [80].

Is suggested to overcome these random flow direction attribution problems for one of the descending neighbors, with the probability proportional to the slope [65]. Other flow direction methods [69, 70] have also been suggested in an attempt to solve the limitations of D8. These attribute a fractional flow to each smaller neighbor, proportional to the slope (or, in the case of the Freeman method, inclination to an exponent) for that neighbor. The multiple flow direction method, here designated by MS (based on multiple slope directions), have the disadvantage that the pixel flow is dispersed for all of the neighboring pixels with lower elevation.

An algorithm was developed using the aspect associated to each pixel to specify the flow directions [71]. The flow is directed as if it was a ball rolling on a plane, liberated from the center of each cell grid. This plane is suited to the elevation of the corners of the pixel, and such corner elevations are estimated by the average of the elevations of the central elevation of the adjacent pixels. This procedure has the advantage of continually specifying the flow direction (an angle between 0 and 2π) without dispersion. Extending the ideas of the previous methodology, a group of elaborated procedures was presented called Demon [66]. Gridded elevation values are used as pixel corners, instead of limiting to the center, and a plane surface is formed for each pixel. The authors recognize the flow as two uniform dimensional origins along the area of the pixel, instead of flow paths drawn from the center of each pixel. The upstream area is evaluated through the construction of detailed flow tubes. It is presupposed that a local plane adjustment for each pixel requires approximation because only three points are necessary to determine a plane. The best adjustment plane, in general, cannot cross the four elevations in the corners, leading to a surface representation discontinuity on the edges of the pixel. The local plane adjustment for specific combinations can lead to inconsistent or deduced flow directions that are a problem in the Lea and Demon methods [63].

Being such, a new procedure to represent the flow directions and calculation of the upstream areas using a grid based DEM was suggested, see reference [63], called infinite D or D∞. The D∞ method calculates the water flow direction according to the steepness of the terrain, distributing the flow proportionally among the neighboring cells. This procedure continually specifies the flow direction (an angle between 0 and 2π) taken in the most accentuated hillside slope, distributing it among the eight facets generate by a 3 x 3 pixel mesh that contains the analyzed pixel in the center. These facets avoid the approximation involved in the plane adjustment and the influence of neighbors with higher altitudes on the upstream water flow [82]. When the direction does not follow one of the cardinal (0, π/2, π, 3π/2) or diagonal (π/4, 3π/4, 5π/4, 7π/4) directions, the accumulated flow is calculated from the flow contribution of a pixel between the two upstream pixels according to the proximity of the flow angle in relation to a right angle for the central pixel. A great advantage of the method D∞ is in considering the form of the divergent surface, in other words, the flow also becomes divergent [83]. Comparing results of the statistical tests and map influence and dependence analysis on the calculation of the upstream area in DEM, a better performance of the D∞ method was observed in relation to the D8, MS and Lea methods, being comparable to Demon, but overcoming its problems of frequent inconsistencies [63].

On defining the drainage directions, it is expected that the resulting drainage network is located within the river channel. Depending on the method used, significant DEM differences in the distribution of the contribution area are obtained. Figure 3 presents the accumulated flow map using the D8 and D∞ methods. As the D8 method routes the whole flow to the cell of higher gradient, rectilinear drainage lines are observed. In turn, the D∞ method, since it considers a proportional distribution among the pixels according to the steepness, does not present the characteristic angular tracings of the flow path restriction.

Using low and high resolution DEM assay data examples, differences are more notable with the increase of the DEM data resolution, especially on the slopes scale [63]. Furthermore, the drainage direction determination methods can produce different results that do not always agree with the reality, mainly when applied in plane areas [60, 84, 85], because they depend on the treatment that each algorithm gives to these regions.

Figure 3.

The accumulated flow by the D8 and D-infinity (D∞) methods. (Source: Elaborated by the authors).

Seeking to define a hydrologically consistent digital elevation model and to obtain a sub-basin drainage network, a study [22] tested the D8 [64] and D∞ [63] methods. In the analysis of the mean error between the observed and estimated drainage, the D∞ method provided higher drainage pathway detailing and agreement among the drainage networks. The D8 method provoked errors in the orientation of the drainage network matrix.

The inability of the single direction method (D8) to simulate the flow direction along the inclination of the hillside was also noted [60], which emphasized best performance of the multiple direction method (D∞). The same observations were made in erosive process spatial distribution analysis studies [20, 83, 86-88]. These results were even verified in a study whose purpose was the modeling of the topographic factor [46], which opted for the multiple flow direction method due to its better adjustment in the erosive process analyses.

3.2.2. Desmet & Govers algorithm

The great disadvantage observed in the USLE/RUSLE models is the two-dimensional evaluation used to determine the effects of the topography. In these models, the landscape has been generically treated as homogeneous, with plane characteristics. The first research that developed a procedure for the soil loss calculation capable to consider the slope form divided the irregular slopes into a limited number of uniform segments was [89]. Continuing this study, weights were attributed for the slope stretches according to their convexity or concavity [3].

Extending the study of [89], the upstream contribution area concept was introduced for the calculation of the L factor [46] that was applied to the RUSLE LS factor equations [90, 91]. For the calculation of the contribution area, a multiple flow direction algorithm was used [70]. The L factor is expressed according to the equation:

$$L_{i,j}=\frac{\left[{\left(A_{i,j-in}+D^2\right)}^{m+1}-{\left(A_{i,j-in}\right)}^{m+1}\right]}{\left[D^{m+2}{x}_{i,j}{}^m({22.13}^m)\right]}$$E8

Where L_i,jis the slope length factor of a cell with coordinates (i, j); A_i,j-inis the contribution area of a cell with coordinates (i, j) (m²); D= is the cell grid size (m); x_i,jis the flow direction value, obtaining the equation x = senα+ cosα, where αis the flow direction angle; mis the coefficient that assumes the values: 0.5, if S≥ 5% (Sthe steepness degree); 0.4, if 3% ≤ S< 5%; 0.3, if 1% ≤ S< 3%; and 0.2, if S< 1%.

For the steepness calculation, the following algorithm was employed [92]:

$$G_{ij}=G_x{}^2+G_y{}^2$$E9

Where G_xand G_yare, respectively, the gradient in the direction x (m/m) and the gradient in the direction y (m/m).

The LS factor for a grid of cells can be thus obtained by the insertion of Li,j and Gij in the LS factor equations of the chosen USLE or RUSLE approach.

This algorithm makes calculations of the steepness, flow direction and the amount of flow that accumulated upstream from a pixel for each pixel [46]. As such, the pixel to pixel topographical factor is calculated along complex slopes. As result, it is possible to define where there is significant distance from the watershed and where there is flow convergence (concave slopes), as well as high steepness, where the LS value tends to be high. In compensation, that value is low in the interfluves (hill tops and plateaus), because the slope length and the steepness are reduced.

The method is not limited to express the sediment transport capacity by the runoff, but it also considers the surface flow, the ramp geometry - if concave or convex - and the erosion way [94]. Comparing the results obtained in the automatic method [46] and manual one [3], a research work [94] verified similar LS values in areas of low steepness. However, in the case of complex slopes, in more sloping areas, the LS values generated by the Desmet & Govers algorithm [46] were significantly superior. That can probably be explained by the assimilation of the convergence and by the respective flow accumulation of this method, which does not occur with the Wischmeier & Smith method [3].

Refining the method of Desmet & Govers, another study [80] incorporated an infinite flow direction model (D∞) and additional methods to isolate the slope length factor. Such method can be applied in situations where a complex topography can influence the surface runoff path and where the excessively long slope lengths, calculated from DEM, can be in need of new landscape detailing. The authors validated the method by the comparison of the statistical distribution of the LS values in GIS with the LS distribution, calculated from field observation data, providing support for applicability of the GIS method to obtain spatial heterogeneity and LS factor magnitude. The evaluation in GIS presented statistical distributions of the LS factor values very similar to those described in the field data, supplying strong support for the use of GIS based methods to represent the spatial heterogeneity and LS factor magnitude for the first time.

3.2.3. RUSLE–3D

From Equation 8, derived by Desmet & Govers, an LS factor equation was generated that is used by the RUSLE 3D model [9]. The model includes irregular hillsides integrating a wide spectrum of hillside convexities and concavities and it incorporates the contribution area Afor the determination of the LS factor:

$$G_{ij}=G_x{}^2+G_y{}^2$$E10

Where A_(r)is the upstream contribution area (m²); β_(r)is the slope inclination angle (degrees); and mand nare flow type dependent parameters.

Typical mvalues are 0.4-0.6 and for n1-1.3. The exponents for the runoff and slope terms in the soil detachment and sediment transport equations reflect the interaction among different flow, detachment and soil transport types. Figure 4 shows the spatial pattern of the topographic potential by RUSLE-3D with different values for the exponent m. For laminar flow (m = 0.1), the detachment and transport of sediments increase relatively little with the amount of water. This type of flow is typical for areas with good plant covering, but also for a severely compacted soil, where the compacting prevents the detachment and formation of furrows. The value of the exponent mof this flow is low, represented by the contribution area [51].

In case the area presents both flow types, usually due to the spatial variability of soil use and properties, an mvalue=0.4 balances the impact of the surface laminar and turbulent flow and it supplies average satisfactory results [51] (Figure 4). When furrow and gully erosion in degraded soils vulnerable to the formation of deep furrows prevails, there are high water flow turbulence conditions and higher water impact, reflecting in a high exponent (m=0.6). The dense vegetation impedes the creation of furrows and maintains a dispersed water flow, while in situations of bare soil, the detachment caused by the flow turbulence increase leads to the formation of furrows [51].

Figure 4.

LS factor of RUSLE-3D with m value of = 0.1 for laminar flow; m = 0.4 for laminar and concentrated flow impact; and m = 0.6 for high flow impact in the erosion pattern. (Source: Elaborated by the authors).

The exponent of the spatial variable based on the covering can vary seeking to increase the negative impact of the disturbed areas and to reduce the impact in vegetated areas [51]. For instance, for forest m=0.2, for pasture m=0.4, for degraded pasture m=0.5 and for degraded areas m=0.6. In a study of these coefficient calibrations for two hydrographic sub-basins with forest vegetation, pastures and native field (prairie) in Mexico [95], the estimated value for m was 0.49. In the evaluation of the water erosion in forest systems [17] the value adopted for the coefficient was 0.4. In a hydrographic sub-basin with prevalence of native forest in Australia, the coefficient m=0.6 was considered more representative [96]. It was observed that the m value is low (m=0.1) when sediment detachment and transport increase relatively little with the amount of water. Thus, the geometric properties of the topography (slope, curvatures) play a more important role in the evolution of the soil detachment and erosion/deposition pattern than the water flow pattern [95].

The RUSLE model supplies an exponent m expressed in function of the slope angle that reflects the predominating dispersed flow in plane hillsides of gentle slope, while the flow in accentuated slopes is more turbulent. However, in the RUSLE-3D model this exponent supplies satisfactory results only for short segments. The formula for m based on the slope can result in values of 0.8 or higher in longer slopes. For slopes with hundreds of meters in length or for the concentrated flow this exponent predicts extremely high erosion rates in RUSLE-3D due to the contribution area [51]. As such, the RUSLE equation for the variable m was developed for the slope length and traditional field applications, and therefore it is not recommended for use as contribution area without evaluating the results through field measurements.

Considering this question, a study conducted the calibration of the m and n parameters through the comparison of the result of each soil loss estimated by the different coefficients, with the soil losses obtained in sample portions inserted in an analyzed sub-basin [48]. Joint analysis the values that presented good performance in relation to the mean error and mean difference for the soil loss estimates were determined. The obtained results, m=0.5 and n=1.0, correspond to those mentioned in [51] (m=0.4 and n=1.0) for areas with high spatial variability of use and soil properties under which both flow types, laminar and concentrated, occur. In fact, the forest use is predominant in the sub-basin and so the laminar flow is favored. At the same time, the high variability of soils and their respective properties also influenced by the relief can generate concentrated flows.

3.2.4. USPED

As the RUSLE-3D model, the Unit Stream Power Erosion and Deposition (USPED) [9] model is derived from USLE and represents its modifications or improvements. The model was developed considering the limitation of the empirical models when estimating the soil loss for convergent and divergent terrain in large areas allied to the Geographic |Information System (GIS). Proposing an adaptation of the contribution area variable the LS topographic factor was derived using the drainage force unit theory to describe the erosive process associated to the laminar and furrow flow in steep hillsides starting from a DEM [55-57].

An advantage of USPED is the fact that it predicts the spatial distribution of the erosion, as well as the deposition rates under conditions of uniform surface flow and high precipitation. Thus, this model can be applied in complex terrains where the erosion is limited by the capacity of the runoff to transport sediment. The topographic index represents the change in the transport capacity of the flow direction, being positive for areas with topographic potential for deposition and negative for areas with erosion potential. The contribution area is used as the representation of the water flow in a place or grid of cells. In USPED, the LS factor equation is [9]:

$$LS=(m+1){\left[\raisebox{1ex}{$A_{(r)}$}\!\left/ \!\raisebox{-1ex}{$22.13$}\right.\right]}^m{\left[\raisebox{1ex}{$\sin {\beta }_{(r)}$}\!\left/ \!\raisebox{-1ex}{$0.09$}\right.\right]}^n$$E11

Where A is the contribution area (m²); θ is the slope angle; and mand nare constants that depend on the flow and soil property types. For situations where the furrow erosion dominates, these parameters are usually established as m = 1.6 and n = 1.3; where the laminar erosion prevails, m = n = 1.0 is considered [57, 97].

In this model, the water erosion in a DEM cell is dependent on the surface runoff in this cell that in turn depends on the upstream drainage area. When substituting the slope length the upstream contribution area generates the erosion network calculated as the convergence of the sediment flow and the deposition network obtained by the alteration in the sediment transport capacity.

Due to USPED computing divergence of the sediment flow, the impact of the exponents is more complex when compared to RUSLE 3D [51]. In USPED the water flow exponent controls the ratio between the erosion extension and deposition, reflecting the fact that the turbulent flow can transport sediments and the impact of the concentrated erosion will be wider than if the flow was dispersed throughout the vegetation. Figure 5 shows the spatial pattern of topographic potential by USPED with different values for the exponent m. For m = 1, a case of dispersed laminar flow and deposition along the hillside. With m = 1.4 we have the case of the influence of both flow types, laminar and in furrows, on the erosion and deposition, with the deposition beginning in the lower third of the hillside and gullying beginning in headwater areas. For m=1.6, furrows and concentrated flows prevail beginning with great force in the headwater areas and turning the erosion even longer and wider with potential for gullying. In this situation, the extension of the deposition areas is even more reduced.

Figure 5.

Spatial pattern of topographic potential for erosion and deposition by USPED with m=1.0; m=1.4; and m=1.6. (Source: Elaborated by the authors).

The variation of the LS factor coefficients of USPED interferes in the soil loss estimates by varying with the relief forms, plant covering and erosive processes. For this reason, the value of the exponents has been documented and established for different climates and areas. For the United States, the value suggested for n in [98] varies from 0.3 to 2.0. In laminar erosion situations the coefficient n=1.0 prevails and where the erosion in furrows is dominant, n=1.3. For a hydrographic basin of forest and agricultural use in Italy, the values adopted were m=n=1 [99], as well as in [95] after the calibration of these parameters. In the identification of the mining impact in an agricultural area of India, the coefficients m=1.6 and n=1.3 were adopted [100], while in Poland, they opted for m=1.4 and n=1.2 for two hydrographic basins with the presence of intense erosive processes [101]. In a sub-basin of forest use located in Brazil, the coefficients m=n=1.0 were determined [48].

As the conceptual models reflect the physical processes that govern the system describing them with empirical relationships, various quantitative evaluation studies of erosion risk in hydrographic basins have opted for the association of the USLE model with an LS factor that reflects the expected surface drainage according to the topography, in order to reach soil loss estimates closer to reality [1, 55, 56, 77, 79, 102-104]. The analysis of the erosion risk in a sub-basin was conducted evaluating the performance of four topographic factor models (USLE, RUSLE, RUSLE 3D and USPED) in the USLE model [48]. The USPED (0.1286 ton ha^-1) and RUSLE 3D (0.0668 ton ha^-1) models did not present statistical differences in relation to the field losses (0.1354 ton ha^-1) and they generated a water erosion distribution meditated by the accumulated flow, while the LS factors of RUSLE (2.74 ton ha^-1) and USLE (3.65 ton ha^-1) overestimated the soil losses [48]. The model considered most efficient in the modeling of the erosion was USPED. This model represented the erosive process in a broad manner when estimating potential erosion and deposition areas, thus allowing to more precisely define the priority areas for conservationist practices under different management sceneries, agreeing with other studies [51, 99]. The advantages of USPED related to the possibility of predicting the spatial distribution of the erosion as well as the deposition rates were stood out in [105], after its comparison with the LS factor of RUSLE 3D. In this context, several research works are opting for the use of the USPED model [9, 45, 51, 98-100, 106, 107].

4. Models in real environmental scenarios

Application these topographic models in the real environmental scenarios have permitted more accuracy and faster estimate of soil erosion in different regions, reliefs and land uses than manual methods. This way, studies have tried to figure out better results applying different LS factor in the water erosion models to estimate and understand this process on watersheds.

LS RUSLE was utilized by [38] in the USLE model to estimate water erosion distribution caused by forest ecosystems in a small watershed and generating soil loss prediction maps according different land use situations. This same LS factor was used in USLE model by [108] allowing identification of the water erosion potential in a watershed forested with eucalyptus and by [109] to estimate of the sediment delivery ratio in a watershed upstream from the hydroelectricity plants.

[110] applied USPED to identify the influence of changing land use on erosion and sedimentation in different land use situations in watershed. [48] applied USLE, RUSLE 3D and USPED in a small watershed for predicting water erosion by eucalyptus plantation founding best results for USPED model. USPED and RUSLE were also applied for assessing the impact of soil erosion/deposition on the archaeological surface at the archaeological site in Greece and USPED presented better results [111]. Using USPED, [112] presented a modeling approach to implement the support practices factor using geographic systems information where data are unavailable and they concluded that USPED with adequate support practices permit to reduce erosion process.

Thus, considering the erosion problems in the world and data available efforts have been done to improve erosion models.

5. Final considerations

The analysis and obtaining of the topographic factors conducted in the digital environment has become a fundamental piece in erosion model progress, because they address the systematic analyses from specific Geographical Information Systems (GIS’s) tools, as well as allow the empirical processing of the data through adaptations of analogical techniques, thus maintaining researcher interpretation. The analysis of the topography in GI enables the analyses of the landscape on a large scale, considers the effects of the topographical complexity more fully in the soil erosion, makes the data processing easier and faster and reduces the relative cost. It is stood out that the reliability of these estimates is directly related to the precision of the topographic surveys used for the derivation of the digital elevation model (DEM).

Evolution of LS factor driven by advanced technology allows the application of different topographic models in the USLE/RUSLE equations modernizing, and improving the estimate of those models. In semi-empirical algorithms, where the contribution area constitutes the central concept, advantages include application in slopes of complex geometries, representation of the surface runoff paths and incorporation of the convergent and divergent flow impact by calibration of empirical coefficients that allow to indicate the qualitative and quantitative effects of the changes in the land use without demanding large spatial and temporal databases. Such models, besides determining the erosion areas on a hydrographic basin level, in the case of the USPED model, even allows to determine the deposition areas, including the erosive process to its full extent.