Predicting soil erosion using Rusle and NDVI time series from TM Landsat 5

– The objective of this work was to evaluate the seasonal variation of soil cover and rainfall erosivity, and their influences on the revised universal soil loss equation (Rusle), in order to estimate watershed soil losses in a temporal scale. Twenty‑two TM Landsat 5 images from 1986 to 2009 were used to estimate soil use and management factor (C factor). A corresponding rainfall erosivity factor (R factor) was considered for each image, and the other factors were obtained using the standard Rusle method. Estimated soil losses were grouped into classes and ranged from 0.13 Mg ha ‑1 on May 24, 2009 (dry season) to 62.0 Mg ha ‑1 on March 11, 2007 (rainy season). In these dates, maximum losses in the watershed were 2.2 and 781.5 Mg ha ‑1 , respectively. Mean annual soil loss in the watershed was 109.5 Mg ha ‑1 , but the central area, with a loss of nearly 300.0 Mg ha ‑1 , was characterized as a site of high water‑erosion risk. The use of C factor obtained from remote sensing data, associated to corresponding R factor, was fundamental to evaluate the soil erosion estimated by the Rusle in different seasons, unlike of other studies which keep these factors constant throughout time.


Introduction
Water erosion is currently one of the main environmental problems for degrading soil and water resources. In addition, it poses a risk to food safety and represents a serious barrier to sustainable development (Telles et al., 2011). Studies on soil use and management in watersheds, especially predictive models, are fundamental to reduce erosive processes.
Models developed to calculate soil erosion can be divided into empirical, conceptual and those based on physical processes (Aksoy & Kavvas, 2005;Kinnell, 2010).
Empirical models such as the universal soil loss equation (Usle) (Wischmeier & Smith, 1978) and its revised version, Rusle (Renard et al., 1997), are used worldwide to provide useful information to support soil and water conservation plans (Kinnell, 2010;Oliveira et al., 2011). Annual soil loss is estimated from the combination of the six factors contained in Rusle: rainfall erosivity (R); soil erodibility (K); slope length (L); slope steepness (S); soil use and management (C); and support practices (P). According to Risse et al. (1993), topography factors (LS) and the C factor exert the greatest influence on the global model efficiency.
The C factor corresponds to soil loss under specific cropping conditions, in relation to that occurring in bare soil (Wischmeier & Smith, 1978). It ranges from 0 (high plant cover) to 1 (bare soil). The Rusle uses five subfactors to calculate the C factor: residual effect of soil use (soil management); soil cover by plant canopy; soil cover by crop residues; roughness of soil surface; and soil moisture (Renard et al., 1997). The evaluation of each subfactor, however, is difficult because of the many possible combinations, and the time spent with data acquisition and analysis (Schönbrodt et al., 2010). Constant values of the Rusle C factor, produced in earlier studies, are usually used to evaluate soil erosion in watersheds. These values, however, do not accurately represent vegetation variation, particularly in large areas, which can result in mistaken estimates of soil loss (Asis & Omasa, 2007). To avoid this problem, some studies have attempted to estimate the C factor from remote sensing information (Asis & Omasa, 2007;Schönbrodt et al., 2010;Bargiel et al., 2013).
The normalized difference vegetation index (NDVI) is one of the main indices used for vegetation monitoring and assessment, which allows the monitoring of the surface spatial and temporal changes. Therefore, from NDVI values, some methods have been developed to estimate the Rusle C factor (Knijff et al., 1999;Asis & Omasa, 2007;Durigon et al., 2014). These methods are developed from regression models which correlate C factor values, either measured in the field or obtained from other studies, to NDVI values generally obtained from TM sensor on Landsat 5, or by the moderate resolution imaging spectroradiometer (Modis) data (Schönbrodt et al., 2010;Le et al., 2012;Alexandridis et al., 2013;Bargiel et al., 2013). However, under tropical climate conditions, the C factor tends to be higher than that calculated by these methods for the same vegetation cover. Therefore, a new method for calculating the Rusle C factor, based on NDVI rescaling, was proposed by Durigon et al. (2014). This method produces realistic values for Rusle C factor and is recommended to be used in tropical regions.
The objective of this work was to evaluate the seasonal variation of soil cover and rainfall erosivity, and their influences on the revised universal soil loss equation (Rusle), in order to estimate watershed soil losses in a temporal scale.

Materials and Methods
This work was carried out from March 2009 to February 2011, at the Palmares-Ribeirão do Saco watershed, located in the south of the state of Rio de Janeiro, Brazil. The watershed is delimited from 22º22'53"S to 22º30'33"S, and from 43°20'50"W to 43º30'15"W, and its topography is mostly rugged. The climate corresponds to Köppen's classification of Cw, with temperatures ranging from 12 to 30ºC, and annual rainfall from 1,100 to 1,700 mm. The area was selected because of the spatial diversity of plant cover and the variety of soil classes and slope.
Orthometric altitude digital elevation model (DEM) of the watershed was obtained from two topographic maps, at a scale of 1:50,000 produced by the Instituto Brasileiro de Geografia e Estatística (IBGE) in universal transverse mercator (UTM) projection. A soil map produced by Embrapa Solos (Centro Nacional de Pesquisa de Solos), based on a semi-detailed survey at a 1:20,000 scale was imported to the ArcGIS, where each mapping unit was reclassified to the second categorical level (suborder), according to the Brazilian soil classification system (SiBCS).
The Rusle model (Renard et al., 1997) was applied to predict soil losses. It uses different procedures from Usle (Wischmeier & Smith, 1978), to determine some parameters: erosivity (R) is calculated in fifteen-day periods; and slope length (L) should be discretized to better characterize steep slope, as follows: A = R K L S C P, in which: A is the annual mean soil loss (Mg ha -1 per year); R is the rainfall erosivity factor (MJ mm ha -1 h -1 ); K is the soil erodibility factor (Mg ha h MJ -1 mm -1 ); L is the slope length factor, dimensionless; S is the slope steepness factor, dimensionless; C is the soil use and management factor, dimensionless; and P means the dimensionless factor for conservation practices.
From 1986 to 2000, the erosivity index (EI 30 ) was obtained using regression equations adjusted to monthly mean rainfall (p) as an independent variable. The EI 30 was calculated for fifteen-day periods by Pesq. agropec. bras., Brasília, v.49, n.3, p.215-224, mar. 2014 DOI: 10.1590/S0100-204X2014000300008 using rainfall data from a meteorological station located in the study area. From 2001 to 2009, data from a pluviograph installed in another meteorological station located in the watershed was used to obtain the EI 30 . In this case, rainfall records showed time intervals of 10 min. Finally, the R factor was determined as the summation of EI 30 for the corresponding time period. Chuveros's software, according to Montebeller et al. (2007), was used to compute the kinetic energy of each rainfall segment using the equation suggested by Foster et al. (1981), and converting data to the international system of units (SI), as follows: Ke = 0.0119 + 0.0873 logI, in which: Ke is the kinetic energy per mm of rainfall (MJ ha -1 mm -1 ); and I is the rainfall intensity (mm h -1 ).
Ke was then multiplied by the rainfall depth in each segment. The values for a same rainfall event were then summed, resulting in the total kinetic energy (Ke, MJ ha -1 ) of each rainfall event. The erosivity index (EI 30 , MJ mm ha -1 h -1 ), calculated for each erosive rainfall event, was the product of the total kinetic energy and maximum rainfall intensity in 30 min (I 30 , mm h -1 ). The sum of many EI 30 indices, previously calculated for each rainfall event, allowed determining erosivity indices for every fifteen-day period. To apply the soil loss model, the adopted R factor corresponded to EI 30 indices calculated for four-month periods, including two months before and two months after image acquisition. This procedure was employed to estimate soil loss in the analyzed 22 images, since soil cover did not exhibit a wide variation over the considered period.
Soil erodibility (K factor) in the watershed, expressed in Mg ha h MJ -1 mm -1 , was estimated according to Handbook 703 of Rusle (Renard et al., 1997). It was calculated for each soil class in the watershed as K = [2.1×10 4 (12 -OM)M 1.14 + 3.25(s -2) + 2.5(p -3)]/ 7.59 × 100, in which: OM is the organic matter content (dag kg -1 ); M is the dimensionless parameter for soil texture; s is the dimensionless soil structure class; p is the dimensionless parameter which depends on the profile permeability (Pe).
Data on soil organic matter, texture and structure were obtained from a semi-detailed survey of soil types conducted in the study area by Embrapa. In addition, soil samples were collected in known geographical coordinates, in order to determine macroporosity using the volumetric ring method.
The M parameter was calculated using equation M = (%silt + %fine sand)(100 -%clay) (Wischmeier & Smith, 1978), and profile permeability (Pe) with equation Pe = Po 2 (Pizarro, 1978), in which: Po is the macroporosity percentage. The p parameter of K factor equation, for each soil profile, was calculated after obtained data in profile permeability equation and according to Wischmeier & Smith (1978), using the permeability classification established by the soil survey division staff (1993).
Two topographic maps, with the 1:50,000 scale of Miguel Pereira (SF-23-Z-B-I-3) and Vassouras (SF-23-Z-A-III-4), were used in this study and included hydrological maps generated by IBGE.
All map layers were transformed from Microstation "dgn" to the ArCGIS "shp" format and imported into the ArcGIS 9.3 software. Then, the units were transformed from kilometers to meters, the layers were georeferenced, and the adjacent vectors of the two topographic maps were joined. All layers of interest, such as hydrography, contours, and land cover, obtained from the various available satellite images, were integrated in the software. Subsequently, the watershed boundaries were demarcated by identifying the drainage divide.
To obtain slope length (L) and slope steepness (S), 10 m spatial resolution DEM, which was previously analyzed to eliminate spurious data, was used. This resolution was chosen by accurately defining flow direction in the watershed, using the "flowdirection" tool of the ArcGIS software. Slope length corresponded to 10.0 m, when water flow was perpendicular to the pixel line, and to 14.14 m when it was diagonally oriented.
The soil use and management factor (C) was determined by using 22 Landsat 5 satellite images (TM sensor), path 217, row 076, from 1986 to 2009. Atmospheric effects on vegetation indices were avoided, by calculating surface reflectance using the 6S model (Antunes et al., 2012). Watershed average altitude (780 m), a tropical gaseous atmosphere, a continental aerosol model, and horizontal visibility (km) of each image were adopted for this atmospheric correction. The images were georeferenced with ground control points, and the NDVI was determined for all images. The NDVI was then used to obtain new images of a rescaled C factor (C r ), as per the following equation (Durigon et al., 2014): Considering the nonexistence of conservation practices in the watershed, the unit value of 1.0 was attributed to the P factor.
Soil loss determined through the Rusle model (Renard et al., 1997), using the ArcGIS software, was classified into four groups of soil loss: 0 to >450.0 Mg ha -1 (Group 1); 0 to >90.0 Mg ha -1 (Group 2); 0 to >9.0 Mg ha -1 (Group 3); 0 to >1.8 Mg ha -1 (Group 4). Mean soil loss was also estimated using mean values of erosivity and soil cover factor for fifteen-day periods. The latter was obtained for 24 fifteen-day periods of each year, using the NDVI from the TM Landsat 5 images above described (Durigon et al., 2014).

Results and Discussion
The differences between the erosivity indices were significant, with the greatest values obtained in images taken in the rainy season, from December to April (Table 1). In Brazil high values of annual precipitation do not necessarily produce higher values of erosivity because of variation in rainfall intensity (Oliveira et al., 2013b). This way, the greatest erosivity values are generally caused by intense rainfall occurring mainly in the rainy season. Further, Carvalho et al. (2012) found similar results of rainfall erosivity for the same study area.
K factor values in the watershed ranged from 0.0120 to 0.0303 Mg ha h MJ -1 mm -1 ( Figure 1A). Soil classes 2 (0.0120 < K < 0.0151) and 3 (0.0151 < K < 0.0154) predominate in the watershed, corresponding to 28.5 and 18.0% of the area, respectively. The highest values were observed in soil classes associated to Ultisol (Argissolo according to SiBCS), except for the Ultisol/Oxisol unit (Argissolo Vermelho-Amarelo/ Latossolo Vermelho-Amarelo, SiBCS). The lowest erodibility values occurred in soil types combining Oxisols (Latossolo Vermelho and Latossolo Vermelho-Amarelo, SiBCS). Ultisols usually show high erodibility values, in contrast to Oxisols. The textural gradient and abrupt texture changes of Ultisols favor erosion because the subsurface horizon exhibits low permeability, facilitating surface runoff of pluvial waters, and ultimately promoting erosion. The more clayey the B horizon, the more difficult the water infiltration and the more severe the erosion process. In an alluvial soil, which consists of complex units of mineral soils, Jebari et al. (2012) recorded erodibility from 0.036 to 0.054 Mg ha h MJ -1 mm -1 .
The topography factor (LS) ranged from 0 to 10.83 (Figure 1 b). Values above 6.0 in 48.5% of the watershed area were found. Areas with the highest LS factor values exhibited the highest slope steepness (44.5 to 146.5%). Some studies have concluded that slope steepness rather than slope length affects the value of the topography factor (Van Remortel et al., Soil cover classes of higher relevance (0 < C r < 0.6), for each image used in the study, were obtained (Table 2). More than 90.0% of the watershed shows C r values less than or equal to 0.3, confirming the predominance of agricultural land use, with permanent plant cover and buildings and bare soil occupying small patches. The C r factor values contrast with the C values of 0 to 0.8767 (Mhangara et al., 2011) and of 0 to 0.5 (Ranzi et al., 2012), which were used to estimate soil erosion in the southernmost portion of South Africa and Northern Vietnam, respectively. This difference is attributed to our approach to determine the C factor (C r ).
The images acquired on 1 October, 1994, and on 2 August, 2007, characterize the dry season in the area, producing the lowest areas with C r values in the first two classes of this factor (0-0.2). The images that revealed the highest areas with plant cover factor in these classes were those obtained on 20 May, 1986, and on 10 May, 1994, both corresponding to the end of the rainy season and onset of the dry season. The highest mean C r values were found for images exhibiting the smallest areas in classes 0-0.1 and 0.1-0.2 (Table 2). These values were 0.235 for images taken on 1 October, 1994, and 0.234 for images from 2 August, 2007. On the other hand, a mean C r value of 0.090 was obtained for images from 20 May, 1986, and 10 May, 1994. It is highlighted that values of 1.0 can be identified only in specific sites of the watershed and do not reflect a relevant spatial condition, mainly on water bodies where the near infrared reflectance is nearly zero, yielding negative values of NDVI.
Despite the similarity between the image groups in relation to soil cover, they differ with respect to relative soil loss, primarily because of rainfall erosivity in each period. The rainfall erosivity value (Table 1), used to estimate soil loss from Figure 2 A, was 1,385.67 MJ mm ha-1 h -1 , and from Figure 2 B, it was 9.50 MJ mm ha-1 h -1 . These data produced respectively mean and maximum losses of 35.46 and 178.50 Mg ha -1 and of 0.23 and 3.10 Mg ha -1 (Table 3). Therefore, the difference between erosivity values produced a substantial variation in mean soil loss. However, soil loss differences between Figures 2 C and D are not as pronounced because their erosivity values are similar in comparison to those of Figures 2 A and B. For 20 May, 1986, and10 May, 1994, mean soil loss was 8.14 and 30.71 Mg ha -1 , and maximum loss was 164.17 and 614.53 Mg ha -1 , respectively (Table 3).
Although Figures 2 A and C correspond to Group 2 (0 to >90 Mg ha -1 soil loss), the greatest losses have been found in the image taken on 1 October, 1994, which exhibits the largest area within the classes of higher soil loss. In these situations, rainfall erosivity and soil cover showed a combined effect on the final outcome, since both were higher in Figure 2 A, which resulted in higher losses. This shows that the initial period of rainy season is critical for erosion protection, as soil is more vulnerable, due to low vegetation density.
Images from Group 1 (0 to > 450 Mg ha -1 ) were obtained in the rainy season (December to May), the period of highest erosivity. In these months, soils usually exhibit a higher moisture content, promoting greater runoff and, consequently, soil loss (Ranzi et al., 2012). For this group, the lowest mean soil losses are associated to images from 10 May, 1994, 11 April, 1995, and 27 April, 2001, which also show the lowest rainfall erosivity in Group 1 (Table 1). Although the highest erosivity value of Group 1 was obtained in the image from 28 January, 1986 (5,394.04 MJ mm ha -1 h -1 ), mean soil loss in this image was lower than that found in the other images of this group because of the good soil cover promoted by vegetation. On 28 January, 1986, the watershed area occupied by the 0-0.1 class of the C factor was 4,252.68 ha (Table 2), which is much higher than the value of 1,106.82 ha in the image from 11 March, 2007, which represents the highest mean soil loss in the group. Results shown for the other groups in relation to mean soil loss are homogeneous, except for the image from 1 October, 1994 (Group 2), whose soil loss is greater than that found for some Group 1 images. This image, however, exhibits a higher rainfall erosivity (1,385.67 MJ mm ha -1 h -1 ) than others from this group (Table 2). Furthermore, the C r factors calculated for this image represented an area of only 119.70 ha (1.4% watershed) in class 0-0.1, whereas for the other images evaluated in this group, the percent areas in class 0-0.1 are 58.7% (20 May, 1986), 51.0% (26 May, 2000) and 26.0% (15 June, 2007). These data explain the difference in soil loss between the 1 October, 1994 image and the others in this group. For the soil loss class ranging from 0 to >90 Mg ha -1 (Group 2), we used images from the end of the rainy season (May, June) and its onset (October). Rainfall in these months usually shows lower depths and intensity, thereby causing lower erosivity. Vegetation in images representing dry season onset (May and June) still shows an adequate soil cover. The highest soil loss was obtained in images from 1 October, 1994, and the lowest on 20 May, 1986. Group 2 images exhibit more homogeneous vegetation than those of Group 1. Images that reflect soil loss in Group 3 (0 to > 9 Mg ha -1 ) are also from May, June and July, which have the lowest rainfall erosivity (154.90 MJ to 339.57 MJ mm ha -1 h -1 ). In other words, because of the previous rainy season, soil cover in the watershed is usually more abundant in this period, promoting a good cover, and consequently decreasing soil loss.
The soil cover has a crucial role on the soil erosion processes in the studied watershed. Some authors have also reported that soil cover has been taken as a key factor to understand surface runoff and soil erosion in different management systems in agricultural lands  or in forest studies (Geißler et al., 2012).
As found in class 2, the lowest soil losses (from 0 to >1.8 Mg ha -1 ) were also obtained in images from May to August, which correspond to the months following the end of the rainy season (Table 2). These losses ranged from 9.50 to 41.18 MJ mm ha -1 h -1 year -1 .
The greatest values for annual soil loss are concentrated in the central watershed area, which holds the greatest classes of erodibility (K) and topography (LS) factors ( Figure 3). The annual soil loss reached 719.97 Mg ha -1 , with mean loss of 109.45 Mg ha -1 .  These studies were developed in watersheds with some similar features, mainly the agricultural land cover (great C factor), and on rugged relief, which consequently promotes great values of the topographic factor. In this study, the soil erodibility and rainfall erosivity factors were estimated using adapted models, and observed soil and rainfall data, while the soil cover factor was from satellite images. The estimates of these factors contain unknown uncertainties, and such study would benefit from time series of soil erosion and sediment yield, which are scarce and difficult to obtain in most cases (Oliveira et al., 2013b). However, soil loss over a long period of time (and a larger area) might be estimated correctly on the basis of the Rusle methodology because overestimations and underestimations can compensate each other, resulting in a good overall assessment of total soil loss (Gabriels et al., 2003). In addition, in the present study, focus is primarily on the evaluation of the erosivity and soil cover factors, and on seasonality effects in estimating soil erosion. It was also proposed a method to be applied in watersheds, aiming at evaluating soil erosion risk areas and contributing to conservation plans for water and soil resources. This way, for these aims the method proposed can be applied in other regions with similar conditions. Conclusions 1. The mean annual soil loss in the evaluated watershed is 109.45 Mg ha -1 ; it reaches nearly 300 Mg ha -1 in the central portion of the watershed.
2. The use of C factor from remote sensing data, associated to corresponding R factor, is fundamental to evaluate soil erosion estimates by the revised universal soil loss equation (Rusle) in different seasons.
3. It is possible to monitor land cover changes, and to evaluate the effects of seasonality on soil erosion estimates through time series of remotely sensed vegetation index.
4. Over the rainy season, rainfall becomes more intense and, thus, more erosive, but soil cover is gradually recomposed with an increase of canopy density, thereby reducing erosion.