A gis-based tool for estimating soil loss in agricultural river basins

Soil erosion is a major problem observed in terrestrial ecosystems. Monitoring and identifying potential areas for erosion becomes extremely important for the better management of these areas. The main aim of this study was to develop a Geographic Information System script tool based on the universal soil loss equation (USLE), which calculates soil loss in three large agricultural sub-basins. Algorithms were implemented in the graphical interface of ModelBuilder and later in Python programming language, thus allowing the creation of a specific script to calculate soil loss in an automatic way. The "USLE Paracatu Watershed" script was validated and proved to be effective in estimating erosion in the three sub-basins with an average processing time of half second per km 2 . This can be added via "ArcToolbox" toolbox in ArcGIS software, so that the user only has to add the variables of the USLE equation and the software will process the algorithms in an automatic way, generating the final map with the soil loss value (t/ha.year). The friendly interface of the script allows it to be used in any area, only requiring the user to enter the updated data of parameters that compose the equation.


Introduction
Erosion processes, terrestrial relief modelers, are natural events taking place worldwide over geological time. However, the propensity of a particular area to soil loss can be enhanced due to its use and occupation, particularly in cases involving the removal of native vegetation (Foucher et al., 2014). Currently, soil erosion is a serious problem that threatens watersheds in tropical regions with intense agricultural use (Sreeja et al., 2015). Such erosion processes not only promote loss of fertile soil, but also increase the amount of particulate matter carried to rivers, reservoirs and seas, with consequent deposition and silting, deteriorating the quality of ecosystems. There is an increasing need for tools that enable efficient and optimized management of these resources. There is no single method that provides assessment and quantification of erosion processes, since erosion calculation is directly linked to spatial scale and the study purpose (Warren, 2002). One of the models most widely used for the calculation of soil loss in agricultural areas is USLE (Universal Soil Loss Equation), developed by Wischmeier (1978) linked to the United States Agriculture Research Service. This seemingly simple equation has shown great potential in the evaluation of erosion in regional scale, both in the identification of potentially erodible areas, as in monitoring erosion over time. Its efficiency is based on the estimation of the main factors that control erosion processes (Wischmeier and Smith, 1978): rainfall erosivity (R), soil erodibility (K), topographical factors (LS), soil use and occupation (C) and conservation practices (P). These factors, when specialized, may be entered into a GIS database, so that the results (output) obtained by geoprocessing techniques are represented in a grid of geo-referenced cells, allowing a quantitative analysis of the spatial distribution of erosion. Soil loss estimated through GIS has lower costs as compared with traditional methods, and greater accuracy for large areas (Lu et al., 2004 andErdongan et al., 2007). Therefore, the use of USLE based on Geographical Information Systems (GIS) is increasingly common, especially in watersheds of developing countries, where agricultural areas are large and resources used for the sustainable management of these regions are scarce (Fistikoglu and Harmancioglu, 2002;Erdogan et al. 2007;Pandey et al., 2007;Dabral et al., 2008;Ozcan et al., 2008;Beskow et al., 2009;Devatha et al., 2015;Guo et al., 2015). Nevertheless, there are few projects that have made use of a tool (script tool) coupled to GIS to optimize this process. In this context, the aim REM: Int. Eng. J., Ouro Preto, 69(4), 417-424, oct. dec. | 2016 of this study was to develop and test a tool on a GIS platform that calculates the erosion rate based on the USLE empirical model. This is a script created through ModelBuilder and ArcPy features (ESRI, 2015). Initially, the conceptual model of USLE was built using ModelBuilder, where spatial problems and techniques to be used were defined, which were used to define input data and structure algorithms. Subsequently, IDE PyScripter was used in the implementation, editing and debugging of algorithms used in the construction of the tool (script tool). The script created was then added as a tool (script tool) to the ArcGIS desktop software 10.2. The developed tool will increase the efficiency in the management of watersheds, reducing the time and allowing its application in areas of different sizes, geomorphologies and uses. In addition to these advantages, being an "add in" in the ArcGIS 10.2 ® software, it becomes a tool with a friendly interface to users, who will not need to be experts in GIS, since they only need to select the input of layers, which in this case are the parameters that compose the USLE formula.

Description of the study area
The Paracatu Watershed (BHRP) is inserted between geographical coordinates 15°30 '/ 19°30'S and 45°10'/ 47°30'W, covering an area of 45,600 km 2 , and 92% of this area, approximately 41,600 km 2 belong to the state of Minas Gerais, and the remaining 5% and 3% are distributed between state of Goias and the Federal District (IGAM, 2006). The climate of the Paracatu Watershed is classified as rainy megathermal Aw type (Köppen. -Pellet et al., 2007), i.e., it presents a summer rainy season from October to April (93% of the annual total rainfall) and a dry winter season from May to October, when rainfall is minimal and the temperature is mild, with an average monthly temperature of 18°C (IGAM, 2006). The study area includes three Paracatu River sub-basins ( Figure 1): the Entre Ribeiros sub-basin (ET), Rio Verde and Rio Escuro. These areas are heavily occupied by mechanized agriculture, in which the irrigation method adopted is by central pivot, in addition to other activities such as livestock and eucalyptus monoculture.

Data sources and processing methods
The satellite images used in this study include SRTM radar imaging (Shuttle Radar Topographic Mission), spatial resolution of 30m, which was used to obtain the LS factor and ASTER satellite imaging to make the map of soil use and occupation of 2006 (C factor). The soil map of the Paracatu Watershed Master Plan (IGAM, 2006), which forms the basis for obtaining the K factor together with historical data from 35 weather stations were used for the generation of the R factor (Figure 1). In a pre-processing step, data were converted to a common projection system (World Geodetic System 1984 Zone 23S). In its original version, the Universal Soil Loss Equation (USLE) is as follows (Wischmeier and Smith, 1978): (1) Equation (2) Equation (3) Equation (4) Equation (5) Equation (6) Equation (7) Where: A is the average annual soil loss (t/ha.year), R is the erosivity factor (MJ.mm /ha h), L is the length factor (dimensionless), S is the slope factor (dimensionless), C is soil use / manage-ment factor and P is the factor related to conservation practices (dimensionless).

USLE factors
The erosivity factor (R) is defined as the sum of the annual average of products of kinetic energy and rainfall by the maxi-mum intensity for a period of 30 min (I 30 ). In this study, the R factor was estimated from historical data of 35 rainfall stations spa-tially distributed over the three sub-basins evaluated in this study (Figure 1), using the equation developed by (Carvalho et al. 1991): El 30 = 111.17 * (r/P) 2 Where: EI 30 corresponds to the average monthly erosivity index (MJ.mm/ha.h), (r) is the average monthly total rainfall (mm), and P is the average total annual rainfall (mm). For determining (R), the monthly values for each rainfall station were added, according to Equation 3: The soil erodibility map (K factor) was obtained by the soil map georeferencing (scale of 1: 500,000) of the Master Plan for Water Resources of the Paracatu Watershed and adapted from (Silva, 2004). The soil use and occupation factor (C) and the conservation practices factor (P) were obtained from ASTER images of July 2006. To obtain the soil use and occupation of sub-basins, a supervised classification method was applied (Likelihood algorithm), using the ENVI 4.2 ® software. This classification resulted in the following classes: 1-Natural vegetation (CP= 0.006); 2 -Eucalyptus plantation (CP=0.012); 3 -Pivots (CP=0.09); 4-Pastures (CP=0.025); 5-Wetlands / water (CP=0.00); 6 -Shade (CP=0.00). The soil use and occupation map was reclassified to the respective C values.
The calculation of the L topographic factor was performed using the algorithm of Desmet and Govers (1996) Where: L ij is the length factor of a cell with coordinates (i, j), A ij-in is the contribution area of the cell with coordinates (i, j) in (m 2 ), D is the cell size (meters), X is the aspect coefficient for the cell grid with coordinates (i, j) and m is the slope function coefficient for the cell grid with coordinates (i, j). Exponent m was calculated from β, which is the ratio between the erosion channels and inter-channels, as suggested by Renard et al. (1997): allowing the automatic calculation of this parameter via ArcGis ® software, resulting in LS maps of the three sub-basins.

Development of the "USLE Paracatu Watershed" script
In the cartographic modeling, the geoprocessing mechanism involves entering input data, performing a mathematical operation in this set of data and returning the results as output data. It was observed that these mathematical operations are carried out with matrix and vector data. The implementation of equations described above (Figure 2) was carried out in the graphical interface of the ModelBuilder in ArcGis ® for desktop software 10.2.2 (ESRI, 2015) and subsequently implemented in Python 2.7 programming language, thereby helping to create a specific script to make the soil loss calculation (A). Then, the code for implementing the Phyton algorithm was created, subsequently used for adding a tool (script tool - Figure 3) in a set of tools (toolbox) in the ArcGIS software with the Spatial Analyst extension. This tool runs in ArcGis versions that support Phyton 2.7. Raster Calculator defines raster cell size based on DEM resolution.
Thus, erosion estimate (A) was obtained from the product (output) of all plans of information related to the USLE factors, using equation 1 (inputs).

Total soil loss (A)
After applying Equation 1, the total soil loss maps (A) were obtained. In order to validate these results, the average soil losses were calculated for each sub-basin, eliminating 10% of their most extreme values, according to the procedure adopted by Silva (2004). Then, the values obtained for soil loss in each sub-basin were transformed into sediment input. The results were then compared to sedimentological data from the Santa Rosa and Porto Alegre stations (Table 1)

Sedimentological data
Unfortunately, there are no sedimentometrical stations in the basins studied, so the erosion values obtained were compared to sedimentological data from two stations of the Paracatu Watershed: Santa Rosa station -covering a drainage area of 13,067? km 2 (downstream of the Rio Escuro sub-basin - Figure 1); and the Porto Alegre station -with drainage surface of 42,171? km 2 (downstream of all subbasins - Figure 1).

Sediment delivery
The total soil loss calculated by USLE is most often greater than the sediment load measured at the mouth of basins (Zhou and Wu, 2008), the re-sult of sedimentation along waterways. The sediment delivery ratio (SDR) is used to correct this reduction. Thus, to validate the results, from values of the average annual soil loss of the sediment delivery in stations, Equation 8 (Lu et al., 2003) and Equation 9 were used (Renfro, 1975).
where SDR is the sediment delivery ratio (%); Y is the sediment supply in the mouth of the basin (t/year); A is the average loss of soil within the basin (t ha -1 /year); a is the contribution area.
The sediment input values in the mouth of each station are in Table 1. In this study, the value was estimated for the year 2006 from the liquid flow x solid flow ratio in the time historical series from 1976 to 2015. The average soil loss calculated for the input area of each station was then compared with value corresponding to that obtained by USLE for each sub-basin, assuming the value obtained as a function of the input area.

Processing time optimization
After installing the script, the processing time is a few minutes per basin, 0.12 seconds per km 2 in computer with processor I7 and 8GB of RAM.

Soil loss and geographical distribution of soil erosions
The three sub-basins show a similar pattern of soil use and occupation. Livestock activity occupies the largest area always above 45%: (49.79% (RV), 52.28% (RE) and 49.45% (ER)), values similar to those of still preserved areas: (43.39% (RV); 41.68% (RE) and 45.89% (ER). The latter situated on the banks of waterways (riparian forests) and slopes (rocky grassland).
Entre Ribeiros: The highest erosion rates were found for Lithic Neosols and Haplic Cambisols (Figure 4), which are associated with cliffs of the São Francisco plateau and Unai ridges, registering annual soil loss values (A) greater than 12.7t/ha.year (Figure 4a).Due to the geomorphology and poor and shallow soils, these areas are considered to be restricted to agriculture. Areas with lower erosion values in the ER basin (Figure 4a) are located in geomorphology unit designated by São Francisco depression, where (A) values range from 0.0t/ ha.year to 12.6t/ha.year. Predominant soils are "red latosols" described as not hydromorphic mineral soils, deep and porous, with thick B horizon, and in addition to these physical characteristics (EMBRAPA, 1999), these are associated with a flat and gently wavy relief, where there are large areas of irrigated agriculture, pastures and reforestation. The lowest erosion values are related to the reduced slope along with the practice of good conservation techniques, despite the intense soil use.
Rio Escuro: Analyzing erosion maps (A) (Figure 4b) it was observed that the highest erosion rates are in Lithic Neosols and Haplic Cambisols, with values above 17.2t/ha.year. Areas of flat and gently wavy relief show the lowest A values, which range from 0.0t/ha.year to 17.1t/ha.year, dominated by red latosols and red yellow latosols, which by their physical and chemical characteristics allow the development of agricultural activities.
Rio Verde: In Rio Verde, the highest erosion rates (A) are found in Lithic Neosols and Yellow-Red Ultisol (Figure  4c), where the relief is more uneven, and thus the LS factor is responsible for the high erosion rates (> 9.4t/ha.year). In the remaining area of the basin, values range from 0.0t/ha.year to 9.3t/ha.year. These low A values are associated to Red Ultisol, Yellow-Red Ultisol and Fluvic Neosol in relatively flat reliefs. These areas are almost exclusively occupied by agricultural activities and pastures. Analyzing the erosion results (A) in the three sub-basins, it was found that for being a region with intense agricultural activity, soil loss in these areas remained within expected values, i.e., between 0.0t/ha.year and 9.3t/ha.year because, according to (Bertoni and Lombardi Neto, 1990), the soil loss tolerance in Brazil ranges from 4.5 to 15.0t/ha.year with average of 10t/ha.year. The highest erosion rates were registered for Lithic Neosols, and these values are related to cliffs. According to soil loss data, although these regions are susceptible to erosion, the soil use and occupation map has shown that these are covered by natural vegetation (forest and savanna), reducing the erosive process. However, if vegetation is removed, erosion can occur.

Results validation
Observing data shown in Tables Table 2 Sediment contribution (%) and soil loss average (A) in the stations and sub-basins.

Conclusions
The results of this study showed that the universal soil loss equation (USLE) was effective in estimating the erosion in large watersheds when adapted to a Geographic Information System (GIS). The "USLE Paracatu Watershed" script developed for these basins allowed the identification and characterization of soil loss in a prompt, effective and robust manner. The processing time for the three sub-basins showed that although areas are large, the average processing time was less than one second per km 2 .
Thus, for future management, this model could be used in a practical and effective way because now the user needs only to enter updated data (DEM, Land Use/ cover map and Rainfall map), and the software will automatically process the algorithms previously described in the methodology, generating the product; in this case, the amount of erosion of the selected area. This tool (script) will contribute to efficiency as it will optimize time and will offer greater flexibility in order to enable the implementation of processes with updated input parameters from different years, or even different locations. According to results obtained by Gomez (2012), the soil erosion modeling on large scales using the universal soil loss equation (USLE) showed good results in the characterization of the Brazilian territory. The automated (GIS) modeling provides an overview of the system, allowing easy access to information by the general public, and therefore the generation of conservation policies in units under study.