Soil Erosion Satellite-Based Estimation in Cropland for Soil Conservation

: Intensive cropland expansion for an increasing population has driven soil degradation worldwide. Modeling how agroecosystems respond to variations in soil attributes, relief and crop management dynamics can guide soil conservation. This research presents a new approach to evaluate soil loss by water erosion in cropland using the RUSLE model and Synthetic Soil Image (spectroscopy technique), which uses time series remotely sensed environmental, agricultural and anthropic variables, in the southeast region of São Paulo State, Brazil. The availability of the open-access satellite images of Tropical Rainfall Measuring Mission (TRMM) and Landsat satellite images provided ten years of rainfall data and 35 years of exposed soil surface. The bare soil surface and agricultural land use were extracted, and the multi-temporal rainfall erosivity was assessed. We predict soil maps’ attributes (texture and organic matter) through innovative soil spectroscopy techniques to assess the soil erodibility and soil loss tolerance. The erosivity, erodibility, and topography obtained by the Earth observations were adopted to estimate soil erosion in four scenarios of sugarcane ( Saccharum spp.) residue coverage (0%, 50%, 75%, and 100%) in five years of the sugarcane cycle: the first year of sugarcane harvest and four subsequent harvesting years from 2013 to 2017. Soil loss tolerance means 4.3 Mg ha − 1 exceeds the minimum rate in 40% of the region, resulting in a total soil loss of ~6 million Mg yr − 1 under total coverage management (7 Mg ha − 1 ). Our findings suggest that sugarcane straw production has not been sufficient to protect the soil loss against water erosion. Thus, straw removal is unfeasible unless alternative conservation practices are adopted, such as minimum soil tillage, contour lines, terracing and other techniques that favor increases in organic matter content and soil flocculating cations. This research also identifies a spatiotemporal erosion-prone area that requests an immediately sustainable land development guide to restore and rehabilitate the vulnerable ecosystem service. The high-resolution spatially distribution method provided can identify soil degradation-prone areas and the cropland expansion frequency. This information may guide farms and the policymakers for a better request of conservation practices according to site-specific management variation.


Introduction
One-third of global land is occupied by livestock (21%) and agriculture (12%) [1]. The increasing demand of a growing population projected to reach 9.6 billion by 2050 from the current 7.7 billion [2] has put the pressure on the Earth's land for food, natural resources, and climate change and created new challenges between humanity and the world's fertile soil resources.
Soil erosion rates caused by the immense expansion of cropland areas have exceeded the rate that soil is formed by one to two orders of magnitude [3]. The process of soil erosion usually involves the detachment, breakdown, transport, redistribution and deposition of sediments [4]. Soil erosion globally has resulted in around 25 to 40 billion tons of increased sediment every year [2]. Land-use change to cropland is responsible for ~80% of the erosion increase. The highest soil erosion occurs in the least developed countries due to the accelerated soil erosion driven by land-use change and poor land management. On the other hand, these countries have the highest potential of soil erosion reduction by conservation agriculture adoption [5].
Sustainable land management (SLM) is treated as an effective land degradation reduction method to achieve the Sustainable Development Goals (SDGs) regarding food, health, water, climate, and land management [6]. SDGs are a global call to end hunger and poverty, protect the planet, and ensure peace and poverty for all by 2030. Despite the SLM's importance for land and soil degradation, identifying large soil management scaling areas in space and time is challenging due to the time-consuming and high cost for data acquisition of the complex climatic, abiotic ecological factors (i.e., soil characteristics and topography), type of land use and land management practices (i.e., tillage and crop rotation). Land-use change does not occur linearly over time, since land rights vary in places and are dependent on the political-economic and legal situation [7].
The Universal Soil Loss Equation (USLE) [8] and its revised version (RUSLE) [9] are the most used empirical models to estimate soil erosion driven by water globally. It aims to guide conservation planning and supports the planners to evaluate and predict the soil erosion rate for each of several alternative combinations of cropping systems and management techniques on any site within the specified limits. The RUSLE equation integrates erosivity, erodibility, topography, cover management, and support practice factors. Furthermore, it considers soil loss tolerance, which ponders the productivity loss due to erosion with the rate of soil formation from parent material, the role of topsoil formation, loss of nutrients and the cost to replace them [9].
Remote Sensing (RS) linked with Geographic Information System (GIS) combined in RUSLE provide data to compute soil erosion with better spatial coverage and accuracy with reasonable costs [10]. The soil spectroscopy advent overcomes the lack of reliable soil data that can be applied to the least developed economies. This method along with the Earth observation-based data can be used to assess the soil erodibility factor based on the environmental soil service (organic matter and mineral composition) with the spectral reflectance data [10,11].
Southeast Brazil is the core of sugarcane production with about 5 million hectares in the 2019/2020 season (62% of the national production) [12]. It was estimated that 600 million Mg yr −1 of annual rainfall soil erosion loss occurs annually in this region [13]. Sugarcane crop in southern Brazil requires a large amount of land to produce biofuels, sugar and electricity. In general, sugarcane production is based on conventional tillage using green mechanized harvesting every year, with a replanting period every five years. The transition from a burned to green harvest system resulted in the maintenance of a large amount of sugarcane straw on the soil surface [14]. The thick layer of straw (ranging from 10 to 20 Mg tons ha −1 ) covering the soil after harvesting resulted in several benefits, such as nutrient cycling [15], carbon storage [16], better soil physical and biological conditions [17,18] and control against erosion processes [19]. However, although straw covers results in several ecosystem benefits, Brazil's more recent sugarcane sector has shown interest in removing part of this crop's residue for bioenergy production, which would lead to increased soil degradation in sugarcane fields [20].
There is a tenuous line between soil degradation and sustainable land management in cropland. In sugarcane fields, straw is an essential source of soil conservation because it benefits soil functioning, e.g., erosion protection cover, soil temperature amplitude reduction, and biological activity increasing [21]. Studies have demonstrated that maintaining 7 to 10 Mg tons ha −1 of straw is enough to cover 100% of the area prevailing the crop's agronomic benefits [22,23]. It suggested that if the straw yield exceeds this amount, it may be used as feedstock for Brazil's bioenergy demand (bioelectricity and cellulosic ethanol).
This paper presents a new approach to evaluate spatial and temporal patterns of soil loss and the impacts of straw removal in the southeastern part of Brazil by integrating environmental variables with the cropland management dataset, field data, RS and GIS data in the RUSLE model. A modern multi-temporal satellite-based estimation method was conducted to assess the soil proprieties and spatial patterns of agricultural land use over the last three decades as well as rainfall data in the last decade. Four soil loss scenarios have been assessed by simulation with four different straw coverage rates of about 0, 50, 75 and 100%, respectively, in five years of the sugarcane cycle (2013-2017) to enlighten the impacts of straw removal changes on soil erosion.
The main novelties of this research are related to the soil erodibility parameters of RUSLE and the soil loss tolerance obtained by the soil attribute maps (texture and organic matter) through the bare soil spectroscopy technique (Synthetic Soil Image). This method combines a Revised Universal Soil Loss Equation (RUSLE) model associated with a time series remotely sensed data to estimate soil loss by water erosion in cropland.

Site Description and Sugarcane Cycle
We evaluate the soil erosion of 500 km 2 of agricultural sugarcane land-use area located in the southeastern part of Brazil, west of São Paulo State (Figure 1a). The erosion processes in this region may deliver sediments to the nearby rivers. Therefore, we expanded soil erosion estimations into 13 sub-basins (~1600 km 2 ) based on the streamflow in order to have a multi-perspective overview of the possible off-site effects of soil erosion on the hydrological ecosystem ( Figure 1a). The region is classified as tropical with dry winter (May-September). The average annual rainfall is 1391 mm yr −1 , and the average annual temperature is 24 °C [24]. The landscape is gently undulating and rolling uplands, with slopes ranging between 0 and 28%. The predominance of sandy-mudstone parental material [25] characterizes the soils' sandy/loam texture dominance. Sheet and interill are preponderating types of erosion in the study site (Figure 1b), but gully erosion sites are common, especially in ascending terrain (Figure 1c).
The sugarcane cycle comprises five years: the plant cane harvest and four subsequent harvesting years (ratoons) (Figure 2). Conventional tillage operations are the most common practices during the plant cane stage and the replanting period, which lies in glyphosate application, subsoiling, plowing, harrowing with physical, chemical, and or biological corrections. Conventional tillage is the most critical period due to the exposed soil period until the crop canopy closure (three to six months later) causing erosion. In areas cultivated with sugarcane, the main attributes of the terrain that contribute to erosion are extreme weather events such as the high intensity and volume of rainfall and soil exposure associated with sloping reliefs. So, the harvest period occurs every year or 18 months, and the maintenance of the straw in the field reduces the risk of soil erosion at this step.

Soil Sampling
Sixty-seven topsoil (Epipedon) samples were collected randomly; 20% of the soil samples validated the soil erodibility maps generated from remote sensing. The soil samples brought to the laboratory were oven-dried for 48 h at 50 °C; then, they were ground and sieved (2 mm mesh) to analyze the soil organic carbon and particle size properties according to [26]. The main soil classes found in the study site, according to the World Reference Base for Soil Resources (WRB) [27] were Gleysols, Ferralsols, Lixisols, Leptosols and Arenosols.

Parametrization of Soil Loss by Water Erosion
Multi-temporal satellite image and field observations assessment provide soil loss estimation. Figure 2 demonstrates the general flowchart of the methodology to obtain soil loss.

SYSI
Cropland use detection digital elevation model used to derive the slope. The Synthetic Soil Image (SYSI) technique based on remote sensing and soil spectroscopy-derived soil attributes maps to calculate soil erodibility, soil loss tolerance and detect the cropland use spatial patterns. The environmental variables (rainfall, erodibility and topography) associated with the cropland patterns and the cropland management dataset (five years of the sugarcane cycle) under different amounts of straw (0%, 50%, 75%, and 100%) estimate soil loss. Soil loss higher than the tolerance indicates high land degradation that contributes to increasing climate change crises, while lower indicates the restoration and rehabilitation of the ecosystem.
The Revised Universal Soil Loss Equation (RUSLE) (Equation (1)) [28] implemented in Geographic Information System (GIS) was used to assess the soil erosion rate (Table 1 and Figure 1). Soil losses by water erosion refer to the amount of sediment that reaches the end of a specified area (cell) on a hillslope [5]. Using a GIS raster scheme means that each cell is independent of the others. Thus, soil erosion is not routed downslope across each cell from hillslopes to the sink area or the riverine systems. The RUSLE model does not capture gullying and other geomorphic processes (i.e., mudflows, landslides and tillage erosion) [29]. However, it indicates the soil loss rate at which gully erosion might be expected to begin [9].
Soil Loss = R * K * LS * C * P (1) The detachment-limited model outputs the mass of soil lost per unit area and time (Mg ha −1 yr −1 ). Sheet, interrill and rill erosion processes are given by the multiplication of six parameters: rainfall erosivity (R-Factor in MJ mm −1 h −1 y −1 ), soil erodibility (K-Factor in tons. h. MJ −1 mm −1 ), slope length (L-Factor in meters), slope steepness (S-Factor in percentage), cropping system (C-Factor), and erosion control practice (P-factor). C and P-factors are dimensionless ( Table 1). The spatial resolutions range from 30 m to 25 km, and the parameters were resampled to a 30 × 30 m cell size to model soil loss as the final output.
The soil loss tolerance threshold (T-value) indicates the maximum rate of soil erosion that can occur and still permits crop productivity to be sustained economically [8,9]. The T value according to [30] is estimated by: (1) soil thickness; (2) soil formation rate; (3) guidelines of the USDA-NRCS; and (4) productivity index. In addition, the T value for a specific soil is a function of the following: the rate of soil formation from parent material; the rate of topsoil formation from subsoil; reduction in crop yield by erosion; soil depth; changes in soil attributes favorable for plant growth; loss of plant nutrients by erosion; the likelihood of rill and gully formation; sediment deposition problems within a field; sediment delivery from the erosion site and; and the availability of feasibility, economically, culturally and socially in addition to soil conservation practices [30]. Cropping management practices with the predicted soil loss rate less than the T-value rate may be projected to deliver less soil degradation ( Figure 2). Thus, sustainable land management practices can reduce the risk of land degradation and mitigate the climate change impacts by promoting the restoration and rehabilitation of the ecosystems, for example, by improving their carbon stock [1].

Rainfall Erosivity Factor (R)
Rainfall erosivity represents the erosive powers of rainfall energy (intensity and duration), a total of rain (volume) and frequency over extended time events. We calculated the R-Factor using ten years (2008-2017) of data of the monthly 3B42 product from the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA). Several studies worldwide have estimated erosivity at fair resolution using the TRMM [11,31,32]. This product provides global coverage of precipitation (mm hr −1 ) with a 3 h temporal resolution and a 0.25 degree spatial resolution [33].
We performed the Modified Fournier Index (MFI) (Equation (2)) using TRMM products to estimate the erosivity factor according to the Erosivity Index (EI30) equations recommended by [34,35]. The EI30 is defined as the product of the maximum rain intensity during a 30 min period. There are three main stations relating EI30 to MFI that fit the study site location (Table 2). Our calculated R-Factor is an average of the three erosivity equations. We used data from seven gauges provided by the Water and Electric Energy Department of São Paulo State (DAEE-SP) located toward the study site to validate the R-Factor obtained from TRMM. P is the average annual rainfall (mm), and Pi is the average rainfall (mm) in the month i.
Here, the topsoil attributes were spatially modeled using an improved methodology with satellite Bare Soil Composite Image [39] based on a Geospatial Soil Sensing System data-mining algorithm. This algorithm flagged bare surfaces from the collection of historical images and aggregated the spatially bare soil fragments into a synthetic soil image (SYSI) [40][41][42]. In this study, we used multi-temporal Landsat satellite images to extract bare soil fragments from 1984 to 2019 and aggregate them in order to obtain all periods of exposed soil for the later creation of a mosaic, which allows viewing the entire area of interest with an exposed soil surface. Landsat imagery has the wavelengths (Bands) of Visible (VIS-Blue/Green/Red), Near Infrared (NIR) and Shortwave Infrared (SWIR1 and SWIR 2) with 30 m of spatial resolution and 16 days of temporal resolution.
MID-Infrared and NDVI indexes were used to remove vegetation and crop residue in each image. We used quality bands to mask clouds and shadows. MID-Infrared represents the normalized difference index calculated from Landsat bands SWIR1 and SWIR2, and NDVI defines the normalized difference vegetation index calculated from Landsat bands NIR and Red. MID-Infrared and NDVI threshold combinations of −0.25 to 0.25 and 0.013 to 0.10. The indexes threshold was determined based on studies related to soil reflectance [40,42] and field observation. We used quality bands to mask clouds and shadows or pixels with inconsistent values, and we choose a filter to use Landsat images with no clouds. The images were aggregated into a single image by the median spectral reflectance value achieving the SYSI of the agricultural spatial patterns from a time interval of 1984 to 2019, with a native spatial (30 m) Landsat product. We used the Google Earth engine platform to acquire the SYSI algorithm.
We performed a partial least square regression (PLSR) method to predict the digital soil mapping attributes by inputting the multi-temporal spectra data from the SYSI and the field measured soil samples (80% employed for calibration and the remaining 20% for validation) at the R program. The models were evaluated based on the coefficient of determination (R 2 ) and the root means square error (RMSE).

Slope Length and Steepness Factor (LS)
The LS-Factor [28] describes the effect of the topography on soil erosion. The L-Factor calculates the slope length, and the S-Factor measures the slope steepness. The L-Factor represents the distance from the overland flow point of origin to the point where the slope gradient decreases, water runoff is streamed into a channel and deposition starts [8,43].
Ref. [44] identify that the S-Factor is not uniform for a whole area; hence, they proposed sub-dividing the slope into several segments. Later, [45] extended this approach to a two-dimensional terrain using the unit-contributing area model (Equations (4)- (6)). Ref. [9] incorporated this approach into RUSLE with slope gradient concepts of [46], and they found that soil loss arises faster in slopes that were steeper than 9% (Equations (7) and (8)).
The LS-Factor was derived from a Digital Elevation Model (DEM) product obtained from the Topodata [44], which filled the original Shuttle Radar Topographic Mission (SRTM) 3 arc-second (90 m) data into an interpolated DEM of 1 arc-second (30 m). We calculated the slope length and slope steepness factors at the ArcGIS platform according to [9,45]. Li,j is the slope length factor for the grid cell with coordinates (i, j). Ai,j is an upslope contributing area for the grid cell with coordinates (i, j) (m 2 ). D is the side length of the grid cell (m). xi,j is a contour length coefficient for the grid cell with coordinates (i, j). m is related to the ratio β of the rill to interill erosion. θ is a gradient of slope in degrees. S is the slope steepness.

Control Practice Factor (P)
The P-Factor measures the ratio of soil loss expected for a specific support conservation practice to the corresponding loss with surface upslope and downslope tillage. The support practices for erosion control usually comprehend contouring, strip-cropping, terracing, and subsurface drainage. These practices influence drainage patterns, runoff velocity, and the direction of water volume concentration [9].
The cropland of the study site mainly applies contour farming in terraces according to the slope gradient. Contour means that farmers implement field practices along contours perpendicular to the normal water flow direction. Consequently, it reduces runoff velocity by increasing the surface roughness, providing more time for infiltration [47], protecting the fertile layer, improving crop yield, increasing income and contributing to the environment. The slope (%) ( Table 3) was obtained by the DEM (30 m) derived P-Factor for the arable cropland. The C-Factor comprises the effect of cropping and management practices on soil loss rates; this factor indicates how the management activities will affect the average annual soil and how the soil erosion potential will be disseminated in time during a conservation plan. The C-Factor value for a particular land use type is the weighted average of Soil Loss Ratios (SRLs) that ranges from 0 to 1 [9]; as it increases to 1, land use degradation stresses high soil threatens. The C-Factor may cause a considerable influence on the erosion calculation, which is defined as the most significant among the RUSLE factors [48].
We used a tool developed for sugarcane crops based on São Paulo State [49], and we calibrated it with the sugarcane stalks yield data from the agroindustry plots from 2013 to 2017. This tool enables many combinations of sugarcane system management as input data to obtain SLR subfactors (Equation (9)) and to model the average annual C-Factor (Equation (10)). It requires geographic location, date of tillage (month), tillage growing (winter, year or 18 months), tillage type (conventional, minimum or no-tillage), covers crop (bare soil, bare fallow, and green crop), crop residue management (0% to 100% of crop residues coverage), number of sugarcane cycles, and management level (low, medium and high).
In this study, we assumed four cover management scenarios to assign soil loss based on the yield average of the five sugarcane cycles and the ratio of 120 kg of dry matter straw per ton of fresh sugarcane stalk [50]. The scenarios are classified in: (i) no straw coverage rate (0 Mg −1 ), (ii) 50% coverage (3.5 Mg −1 ), (iii) 75% coverage (5.25 Mg −1 ), and (iv) 100% coverage (7 Mg −1 ).
SRL is calculated for a given condition using five subfactors and has to be calculated for different periods, as [9] recommended. Each subfactor contains cropping and management variables that affect soil erosion. PLU (Prior-Land Use) expresses the influence on soil erosion of subsurface residual effects from previous crops and the effect of previous crop management practices on soil consolidation. CC (Canopy Cover) indicates the effectiveness of vegetative canopy in reducing the energy of the rainfall striking the soil surface. SC (Surface Cover) affects erosion by reducing the transport capacity of runoff water by causing deposition in ponded areas and decreasing the surface area susceptible to raindrop impact. SR (Surface Roughness) is a function of the surface's random roughness. SM (Soil Moisture) influences infiltration, runoff, and soil erosion [9]. Here, soil moisture reflects the soil field capacity, value 1, due to their intensive condition under the erosive process.

Soil Loss Tolerance
Various methods to assess soil loss tolerance have been recommended by soil scientists worldwide; each approach has its assumptions, advantages, and limitations [51]. Land cover and carbon stocks are the principal indicators to estimate land degradation because they can change rapidly [52]. Here, the spatial modeling of soil loss tolerance was estimated based on [53] (Equation (13)), which considers the soil bulk density method obtained in [54] (Equations (11) and (12)) that depends on SOM and soil drained capacity. The soil class map of the agroindustry determined the soil drained capacity. The SOM map was estimated by the satellite bare soil surface methodology (Item 2.3.2). A soil depth of 0 to 0.30 m, and the temporal constant founded on the concept that forming 1000 mm of soil required 1000 years [55], were incorporated in the calculations to obtain T.
where BD is soil bulk density (g cm −3 ), and SOM is organic matter (%) BD is soil bulk density (kg m −3 ). SOM is soil organic matter (g kg −1 ). T is soil loss tolerance (Mg ha −1 yr −1 ), H is soil depth (m), and 1000 is the temporal constant founded on the concept that forming 1000 mm of soil required 1000 years [53]; a 10,000 factor was applied to transform the data from Mg m −2 to Mg ha −1 .

Rainfall Erosivity Factor
The spatial data average rainfall erosivity calculated represents 6078.00 MJ mm −1 ha −1 yr −1 (Figure 3), which according to [34,54] indicates strong erosivity, 300% higher than the global average [5]. We calculated the R-Factor using TRMM data between 2008 and 2017 and the erosivity equations obtained by the daily rain data as a function of the Modified Fournier Index (MFI) (Figure 3). We compared our estimated R-Factor to that calculated from the monthly rain data gauges covering the sub-basins. In the southeastern part of Brazil, the rainy period marks are that warm temperatures start in September and end in May ( Figure  4a). The relation has an R 2 = 0.90, demonstrating a good correspondence between the values with a close 1:1 relationship (Figure 4b).

Soil Erodibility Factor Obtained from the Digital Soil Attributes Mapping
The summary statistics of the soil property in Table 4 show soil attributes physical and organic matter variability at 0 to 0.30 m depth ( Table 4). The low amount of clay and SOM content and high sand content ranging from 73% to 92% demonstrate soil erosionprone features and the low natural fertility in the soil.  The geology map [25] showed that the sandy-mudstone lithology of the Rio do Peixe  Table 5) [55]. By approaching the rolling uplands, the sandstone from the Marilia formation arises. Basic intrusive diabase rocks occur between the sand-mudstone of the Serra Geral formation neighboring the site, which provides the reddish color to the soil types and clay accumulation in the Lixisol subsurface ( Figure 5). Ferralsols, Lixisol, and Gleysols represent 40%, 58%, and 2% of the land, respectively.
Soil permeability was established based on soil type, color and texture ( Figure 5 and Table 5). The permeability decreases as color changes from reddish to yellowish and from sandy to clay texture. We classified the soil classes, in dry conditions, from the agroindustry map and the regional map [56] at the suborder level according to the Brazilian Soil Classification System (SiBCS) [57] and the corresponding classes of the World Reference Base [58].  The multi-temporal satellite bare soil resulted in 35 years of Landsat data from 1984 to 2019 (Figure 6a). The synthesis of all images composes a final image of bare soil, covering 72% of the cropland surface (Figure 6b). The soil surface corresponds to the agricultural inventory data, in which there is a specific time window to monitor soil surface satellite imaging due to the conventional crop management tillage adopted. The remaining sites correspond to forest or grass assigned as NA values for modeling. Figure 6. (a) Multi-temporal satellite images masks using NDVI and Mid-Infrared, (b) Soil Synthetic Image has shown the agricultural land in false color (RGB 543) of the Landsat data with the soil samples sites location used for the soil map attributes prediction by Partial Square Least Regression, (c) Soil spectral profile average from the soil samples, (d) Soil attributes map (clay, fine sand, silt, and organic matter) used to derive the erodibility factor map. Soil texture analyzes resulted in the following classes: Sand, Coarse Sand, Fine Sand, Silt and Clay. Note that only Clay, Fine Sand and Silt maps were generated to calculate erodibility, so the final percentage value will not be 100%. This would require all maps of all textural classes to reach 100%.
The spatial-spectral patterns are related to the soil mineralogy, granulometry and organic matter content identified throughout the satellite image's color and the spectral data (Figure 6c). Bright colors indicate low soil organic carbon content and higher quartz proportions [40]. The shape and intensity of the spectral profile discriminate sandy soil as it increases from the higher reflectance and the ascending shape of the wavelength from Band 1 to Band 5 [39].
The spectral signatures of the surface reflectance provided from the site-specific soil samples validate the bare soil composite image's reliability and quality. Meanwhile, the validation prediction of clay, fine sand, silt, and organic matter maps corresponds to R 2 of 0.67, 0.59, 0.68 and 0.55 and RMSE (g kg −1 ) of 25.8, 48.3, 4.38, and 1.99, respectively, which indicates good agreement and low error between the observed and the predicted dataset.
The spatial distribution of K-Factor values increases where PVA and PV occur in rolling uplands, as observed in the sub-basin 10 ( Figure 7). In these sites, the sand-loamy textures soil particles are easily detached and transported by overland flow. The mean K-Factor value for the agricultural study site corresponds to 0.019 tons h MJ −1 mm −1 . Predicted versus measures K-Factor values present satisfactory performance with an R 2 of ~0.78 (Figure 8).

The Topographic Parameters and Control Practice
The LS-Factor derived from flow accumulation, flow direction, and slope using a Digital Elevation Model (DEM) performed by GIS incorporates surface runoff into soil erosion (Figure 9). The L-Factor stretches the impact of slope length while the S-Factor delivers the effect of slope steepness. The DEM captured the topography changes with precision and estimated soil erosion with accuracy. Ref. [45] demonstrated that the LS-Factor model is suitable for landscape-scale soil erosion modeling. The P-Factor map depended on the slope derived from DEM and contouring as specific support management (Figure 8). The P-Factor result presents low variation; 88% of the area has represented a value from 0.5 to 0.6 with an average of 0.59.

Cover Management Factor
The spatial C-Factor of crop residue management for the scenarios expresses the land use management applied from 2013 to 2017 ( Figure 10). It generated 98 combinations of sugarcane management in the cropland plots. We can observe in sub-basin 7 the cropland plots contrasting two tillage combinations. The sugarcane cycles range from one to nine in the unit plots; the mean is five crop cycles (Figure 11). During the sugarcane crop cycle, plant cane had higher production and yield decay over the next years (Figure 11b). We can observe that about 80% of the area is productive during the year, while the remaining area is under replanting (Figure 11a).  During the sugarcane replanting period, the crop residues are incorporated into the soil by tillage operations, and the soil is exposed, characterizing the prior land use subfactor as bare soil for all the cropland plots. The average C-Factor values from the cropland plots for 0%, 50%, 75%, and 100% coverage are 0.12, 0.0975, 0.0921 and 0.0884, respectively. We applied this C-Factor value for the entire sub-basins dimension.

Soil Loss in Agricultural Regions
The soil loss rates from high-resolution spatially distributed modeling (30 × 30 m cell size) were estimated by integrating the erosivity, erodibility, topography, conservative practices, and cover management data. Our results showed that soil erosion dynamics under different straw coverage rates ( Figure 12) decreases soil loss rates by 22%, 26% and 28% in all sub-basins as the amount of straw coverage increases to 50%, 75%, and 100%, respectively (Table 6).  The average bulk density of 1.69 Mg m −3 derived soil loss tolerance (T-value); the Tvalue presents an area-specific average of 4.3 Mg ha −1 . Our data demonstrate that the annual average soil erosion exceeds the T-value threshold in 12 of the 13 sub-basins under no cover of straw (Table 6). Sub-basins 3, 6, and 7 show soil loss per percentage of the land lower than the T-value (Figure 12a).
The most intensively eroded region were identified in sub-basins 4, 8, 9, 10, 11, 13 with an annual average soil erosion higher than the T-values in sites 100% covered with straw (Table 6). Sub-basin 4 presented the lowest values of total soil loss (52.65 Mg × 10 3 yr −1 ) (Table 4), which was due to the smaller size of the sub-basin compared to the others. The average soil loss scenario of 100% straw coverage is equivalent to the T-value ( Table  6), indicating that all the amount of straw is necessary to balance soil loss tolerance from the perspective of soil erosion. Thus, it is equivalent to 7 Mg −1 of straw according to the sugarcane yield data documented ( Figure 11). However, this amount of straw has not been satisfactory to protect the soil against erosion in 40% of the agricultural land, which was classified with a higher potential of soil erosion estimated at ~6000 Mg 103 y −1 ( Figure  12a, Table 6). We estimated an overall increase in soil erosion amount of 4%, 10%, and 41% under 75%, 50% and no coverage, respectively, which was driven by the possibility of spatial removal.

Discussion
The soil is vulnerable to erosion by nature according to a slow constructive process of fertile soils caused by geology, topography and climate factors. In contrast, anthropogenic land-use changes and unsustainable land management in agriculture can accelerate land degradation [4].
In this paper, we assessed a multi-temporal soil erosion pattern driven by water in sugarcane fields using the RUSLE model at 30 m resolution. We used remote sensing data and GIS techniques to extract R-, K-, and LS-Factors and an exhaustive set of C-Factors, an accurate DEM to compute LS-and P-Factors, and an innovative digital soil mapping method to calculate the K-Factor. Thus, the variability demonstrated by this method has a considerable potential to identify hotspots and concern areas for both uses in the cropland plots scale and the whole area for conservation planning (Figure 12).
Rainfall intensity emphasizes the region's climatic vulnerability, resulting in high values of erosivity factor (Figures 3 and 4), which is a critical climatic driver of soil erosion. According to [31], erosivity correlates with total rainstorm energy to the storm's maximum 30 min of rainfall intensity. The spatial variability of erosion power of rainfall tends to increase or decrease in various combinations driven to the degree of global warming caused by climate change [61], which impacts directly on land degradation, resulting in increasing temperatures, changing rainfall patterns, and intensification of rainfall [1].
Individual powerful rainstorms of a short period can remove topsoil during a rainstorm event, initiating a gully and mudslide, as revealed in Figure 1. Damages caused by heavy rainfall events have a century-scale nutrient depletion, soil particle loss, and highwater treatment costs. The 3-hourly TRMM data provided a good indicator of high-intensity rainfall events (more than 30 mm), which were typically captured at the beginning of the rainy season [62].
Hydrological changes affect soil fertility due to topsoil layer removal, water availability, and the vulnerability of smallholder and subsistence agriculture [63]. Each soil type has properties that have different resistance to the raindrop and runoff. High rainfall erosivity increases soil loss rates that deplete the soil structure (breakdown aggregates) and accelerate the decomposition of the organic carbon matter by a microbial process [4]. Controversially, a decrease in rainfall erosivity may enforce agribusiness development [64] and smallholders, especially in regions with crop irrigation liability. Meanwhile, a decrease in rainfall erosivity in the southeast region of Brazil may suggest a favorable scenario for the continued sugarcane expansion [64] since sugarcane does not need irrigation [65].
The projection of climate change in rainfall erosivity is essential, but at this stage, developing strategies to adapt to these changes is the key to reduce vulnerabilities and improve resilience; e.g., public policies that focus on soil and water conservation such as sustainable agriculture and conservation practices must be encouraged, ensuring food security and energy.
Our digital soil attributes maps obtained extracted from satellite data are instruments for sustainable land use ( Figure 6). Previous studies based on Landsat data have also derived reliable soil proprieties maps demonstrating the substantial potential of these products for the land management applications [66][67][68][69]. We identified that sites with moderate soil loss are mostly homogeneous with topsoil texture that varies from loam to sandy. The loam texture sites demonstrate the occurrence of Ferralsols ( Figure 5) in the flat position of the landscape. On the other hand, high soil loss occurs mostly in Lixisol ( Figure 5) that is over smooth to gently undulating relief with the textural gradient of sandy/loamy at the surface and loamy/clayey in subsurface layers ( Figure 6), featuring the highest values of LS-Factor (Figure 9). We confirmed [70] assumption that an increase in the textural gradient leads to a decrease in soil loss tolerance and an increase in soil K-Factor (Figure 7). SYSI data can capture the agricultural land-use change. Here, we could enlighten two time periods to investigate land-use change using SYSI: (1) from 1984 to 2009, the period with land changes from forestation to sugarcane under harvest burning system adoption; (2) after 2009, with the establishment of the current machine harvest system characterized by the remarkable effect of soil conservation agriculture provided by straw on topsoil. On the other hand, straw removal practices may threaten sugarcane yield, compromising bioenergy production [71]. The nutrient addition in soil in the first period, associated with straw removal of the second period plus the climate change impacts related to land use cover, may be responsible for soil organic matter depletion [16,71]. Our investigations reveal that the high erodibility values in sub-basin 10 correspond to Lixisols, which are associated with the low organic matter content (Figure 7). Lixisols when exposed due to changes in land use and/or soil management lead to erosion processes. These soils have a more clayey B horizon, which reduces water infiltration, favoring surface runoff and, consequently, erosion.
Strategies for soil conservation as cover management are the main component to control soil erosion runoff potential, maintain soil structure and conserve soil organic carbon with low cost [72]. Crop residue retained in the field delivers maintenance to the soil quality; it affects sugarcane yields in all soil types in different magnitudes. Straw can help maintain soil moisture and mitigate water deficit, which stimulates the microclimate by regulating the thermal amplitude and improving the soil biota [73,74]. The yield annual average, inferior to the national average yield of 76 Mg ha −1 [75] (Figure 11a) aggregated with soil compaction, reveals low soil fertility. Field observations and the literature indicated high soil compaction in this region, particularly if straw has been removed [17], which resulted in low soil permeability capacity.
Ferralsol is very well-developed, deep and unsaturated, and it is characterized as having fast permeability. However, this soil was classified as moderate to fast or moderate. We also observed that permeability decreases as the color changes from reddish to yellowish and sandy to clay texture (Table 5). Whereas, as permeability decreases, Lixisol is more vulnerable to topsoil loss, which is well-developed with a textural gradient in the subsurface featuring erosion-prone soils [76]. These attributes result in low yielding potential and straw production. It explains the inefficiency of 100% straw coverage capacity of holding gross erosion and preventing soil degradation in many sub-basins, especially in sub-basins 4 and 8 that presented the soil erosion hot spots (Figure 12d and Table 6) even under a low heterogeneity topography ( Figure 9).
Our C-Factor results demonstrate that the number of sugarcane cycles and the planting season are the SRL that have the most significant impact on the management cover, which is consistent with [49]. The sugarcane cycle can influence the crop residue maintenance in the erosion process by surface cover, canopy cover, and tillage effectiveness in the yielding. While the planting season is affected by the climate, the soil is affected by the following conditions: (i) during replanting, when the soil is entirely uncovered; (ii) at the beginning of the sugarcane ratoon, the previous cycle removed the straw, and the sugarcane canopy cover is not entirely close; (iii) the low number of sugarcane cycles (less than five due to the low yield and coverage potential; and (iv) tillage management operations during the rainy season.
Alternative management strategies on-site for more sustainable sugarcane production are required to compensate for the adverse erosion effects on soil, water, and biotic balance as vegetative barrier conservation practices [77]. Techniques to improve soil organic carbon and increase yield include reduce or no-tillage, crop rotation (minimum soil disturbance) [16], legume cover crop [78], organic fertilizer, filter cake, ashes, biochar, and others. Reduce tillage and no-tillage system concepts were implemented in this century. A reduced tillage system associated with a part of the straw cover in the soil could enhance the soil organic carbon stock, sustain sugarcane yield over the crop cycle, and allow part of the straw to be used for bioenergy [16]. Each agricultural system affects short-term soil CO2 emissions. Undisturbed soil keeps higher soil moisture than conventional tillage; moisture is a control temporal variability factor of CO2 emission. The adoption of no-tillage in sugarcane areas would prevent 30% of soil CO2 emission in tropical soils compared with conventional tillage (La Scala et al., 2006).
The erosion cost could be reduced by 81.2% adopting no-tillage, while the production costs increased only by 0.47 for the soybean crop [79]. Another study observed the same crop and an erosion cost reduction from the different management systems, including notillage, reduced tillage and conventional tillage of USD 15, 16, and 25 ha −1 yr −1 , respectively [80]. Soil erosion assessment is essential from an economic perspective since it degrades soils on-site, producing loss of fertility, and reduces water storage capacity, compromising yield, which over the long term may depreciate land value [81,82]. The hydrographic basins unit is vital to assess the cost of soil loss due to the off-site impacts of land degradation. Soil erosion effects in the surrounding areas can be severe to the freshwater systems causing sedimentation and eutrophication, enhancing urban areas flooding, and impacting the marine ecosystem.
The results reported herein indicate that conservation-effective measures to address on-site erosion are essential to reduce or reverse soil degradation and minimize the climate change impacts from off-site erosion. Furthermore, the obligation of soil conservation practices may help to internalize the costs of the land user. Fertile soil is a limited and non-renewable natural resource, and soil degradation impacts the world from many aspects, e.g., socio-economic, social, political, and cultural for current and future generations.

Conclusions
The method to estimate soil loss indicated soil degradation-prone areas that require sustainable land management. The bare soil surface image obtained from multi-temporal satellite images covering 100% of the agricultural land use is equivalent to 72% of the study site. The image's spectral patterns presented an accurate spectral quality that permits capturing the soil interactions, featuring sandy soil and low organic matter content properties. Our spatiotemporal erodibility factor reveals the correlation of soil organic matter depletion with high erodibility values in Lixisol, which resulted in low soil loss tolerance.
The bare soil surface technique indicates that intensive farming caused high soil exposure rates in the last three decades. Erosion intensification has contributed to land degradation for at least 20 years from the beginning of the multi-temporal series until 2010, when bare soil frequency decreased due to the gradual change to agricultural conservation systems established in the sugarcane fields. Our result indicates vulnerable areas to hold gross erosion in locations with 100% crop residue coverage. This susceptible area is a combination of natural erosion-prone areas with low effectiveness of the tillage practices and the number of cycles (inferior to five) that result in low sugarcane straw yield, which is easily and quickly decomposed by the microbial process. In addition, we suggest further studies using SYSI to deepen soil erosion research.
Furthermore, the strong rainfall erosivity of the region reinforces soil erosion by water, especially in the planting season. Public policies that focus on soil and water conservation to ensure food and energy security are strategies to reduce the vulnerabilities and improve the resilience of the environment to adapt to the increasing temperatures, the changing of the rainfall patterns, and the rainfall intensification driven by climate change. This technique provides multidisciplinary uses that extend beyond sustainable agriculture developments and land-use change monitoring as the costs to replace nutrients loss derived from soil erosion and carbon stock dynamics are the foundation for the emission levels analyses.
The high-resolution spatially distribution method provided can identify soil degradation-prone areas and the cropland expansion frequency. This information may guide farms and the policymakers for a better request of conservation practices according to sitespecific management variation.