Improving the global applicability of the RUSLE model – adjustment of the topographical and rainfall erosivity factors

Large uncertainties exist in estimated rates and the extent of soil erosion by surface runoff on a global scale. This limits our understanding of the global impact that soil erosion might have on agriculture and climate. The Revised Universal Soil Loss Equation (RUSLE) model is, due to its simple structure and empirical basis, a frequently used tool in estimating average annual soil erosion rates at regional to global scales. However, large spatial-scale applications often rely on coarse data input, which is not compatible with the local scale on which the model is parameterized. Our study aims at providing the first steps in improving the global applicability of the RUSLE model in order to derive more accurate global soil erosion rates. We adjusted the topographical and rainfall erosivity factors of the RUSLE model and compared the resulting erosion rates to extensive empirical databases from the USA and Europe. By scaling the slope according to the fractal method to adjust the topographical factor, we managed to improve the topographical detail in a coarse resolution global digital elevation model. Applying the linear multiple regression method to adjust rainfall erosivity for various climate zones resulted in values that compared well to high resolution erosivity data for different regions. However, this method needs to be extended to tropical climates, for which erosivity is biased due to the lack of high resolution erosivity data. After applying the adjusted and the unadjusted versions of the RUSLE model on a global scale we find that the adjusted version shows a global higher mean erosion rate and more variability in the erosion rates. Comparison to empirical data sets of the USA and Europe shows that the adjusted RUSLE model is able to decrease the very high erosion rates in hilly regions that are observed in the unadjusted RUSLE model results. Although there are still some regional differences with the empirical databases, the results indicate that the methods used here seem to be a promising tool in improving the applicability of the RUSLE model at coarse resolution on a global scale.


Introduction
For the last centuries to millennia soil erosion by surface runoff has been accelerated globally due to human activities such as deforestation and agricultural practices (Bork and Lang, 2003).Accelerated soil erosion is a process that triggers land degradation in the form of nutrient loss, a decrease in the effective root depth, water imbalance in the root zone and finally also productivity reduction (Yang et al., 2003).It is widely recognized that soil erosion has been a major threat to sustainable agriculture and food production across the globe since the start of agricultural activities (UNCCD, 2012;Walling, 2009).These effects of soil erosion are currently exacerbated by the global population growth and climatic changes.Organizations such as the United Nations Convention to Combat Desertification (UNCCD) try to address this problem by stating a new goal for Rio +20 of zero land degradation (UNCCD, 2012).
Another aspect underpinning the relevance of soil erosion on the global scale is the effect of erosion on global nutrient cycles.Recently, the biogeochemical components of Earth system models (ESMs) became increasingly impor-Published by Copernicus Publications on behalf of the European Geosciences Union.
tant in predicting the global future climate (Thornton et al., 2007;Goll et al., 2012).Not only the global carbon cycle but also other nutrient cycles such as the nitrogen and phosphorous cycles cannot be neglected in ESMs anymore (Goll et al., 2012;Gruber and Galloway, 2008;Reich and Hungate, 2006).Soil erosion may have a significant impact on these biogeochemical cycles through lateral fluxes of sediment, but the impact on the global scale is still largely unknown.For example, Quinton et al. (2010) showed that erosion can significantly alter the nutrient and carbon cycling and result in lateral fluxes of nutrients that are similar in magnitude as fluxes induced by fertilizer application and crop removal.Regnier et al. (2013) looked at the effect of human-induced lateral fluxes of carbon from land to ocean and concluded that human perturbations, which include soil erosion, may have enhanced the carbon export from soils to inland waters.
In general, the effect of soil erosion on the global carbon cycle has received considerable attention after the pioneering work of Stallard (1998), who proposed that global soil erosion can result in sequestration of carbon by soils.After his work, the effect of soil erosion on the carbon cycle has been studied extensively, but there remains a large uncertainty in the effect of soil erosion on the carbon cycle.For example, several recent global assessments of the influence of soil erosion on the carbon cycle indicate a large uncertainty with a range from a source of 0.37-1 Pg C year −1 to a net uptake or sink of 0.56-1 Pg C year −1 (Van Oost et al., 2007).Thus, in order to better constrain the global carbon budget and to identify optimal management strategies for land use, it is essential to have accurate estimates of soil erosion and its variability on a global scale.
Currently, there exists a large uncertainty in the global soil erosion rates as can be seen from recent studies that show rates between 20 and 200 Pg year −1 (Doetterl et al., 2012).This indicates that modelling soil erosion on a global scale is still a difficult task due to the very high spatial and temporal variability of soil erosion.Different approaches were previously applied to estimate soil erosion on a large or global scale.Most of these approaches are based on extrapolated data from agricultural plots, sediment yield or extrapolated river sediment estimates (Milliman and Syvitski, 1992;Stallard, 1998;Lal, 2003;Hooke, 2000;Pimentel et al.,1995;Wilkinson and McElroy, 2007).
An alternative approach is based on the use of soil erosion models, in order to be able to predict soil erosion rates for the past and future.One of the most applied models to estimate soil erosion on a large spatial scale is the semiempirical/process-based Revised Universal Soil Loss Equation (RUSLE) model (Renard et al., 1997).This model stems from the original Universal Soil Loss Equation (USLE) model developed by the USDA (US Department of Agriculture), which is based on a large set of experiments on soil loss due to water erosion from agricultural plots in the USA.These experiments covered a large variety of agricultural practices, soil types and climatic conditions, making it a potentially suitable tool on a regional to global scale.The RUSLE model predicts the average annual soil erosion rates by rainfall and is formulated as a product of a rainfall erosivity factor (R), a slope steepness factor (S), a slope length factor (L), a soil erodibility factor (K), a land cover factor (C) and a support practice factor (P ).The RUSLE model was first applied on a global scale by Yang et al. (2003) and Ito (2007) for estimating the global soil erosion potential.Various limitations were observed when applying this model on global scale.Firstly, the model is originally developed to be applicable on the agricultural plot scale.This makes the model incompatible with the coarse spatial scale of global data sets on soil-erosion-influencing factors such as precipitation, elevation, land use and soil characteristics.Secondly, the RUSLE and USLE models were parameterized for environmental conditions of the USA and are thus not directly applicable to other areas in the world.Thirdly, only sheet and rill erosion are considered.Finally, the RUSLE model does not contain sediment deposition and sediment transport terms, which are closely linked to soil erosion.
However, the RUSLE model is to our knowledge one of the few erosion models that has the potential to be applied on a global scale due to its simple structure and empirical basis.Therefore, it is of key importance to address the abovementioned limitations first.
To address the first two limitations, Van Oost et al. ( 2007) presented in their work a modified version of the USLE model for application on agricultural areas on global scale.They based their model on large-scale experimental soil erosion data from the USA (National Resource Inventory, NRI, database; USDA, 2000) and Europe by deriving reference factors for soil erosion on agricultural land and for certain USLE parameters.They also introduced a procedure to scale slope, which is an important parameter in the topographical factors S and L of the USLE/RUSLE model.In this scaling procedure slope was scaled from the GTOPO30 1 km resolution digital elevation model (USGS, 1996) to the coarser resolution of the erosion model.This method was based on high resolution OS (Ordnance Survey; 10 m resolution) and SRTM (Shuttle Radar Topography Mission) data on elevation (90 m resolution, International Centre for Tropical Agriculture, CIAT) for England and Wales.Doetterl et al. (2012) showed that together with the S factor, the rainfall erosivity or R factor explain up to 75 % of the erosion variability across agricultural areas at the large watershed scale.These factors represent the triggers for soil erosion by providing energy for soil to erode.They can also be seen as the natural components of the RUSLE model, as they include very little or no modification by human activities (Angulo-Martínez et al., 2009) apart from indirect effects on precipitation and extreme events due to anthropogenic climate change.In this way they represent the natural environmental constraints to soil erosion that are important to capture before the effect of human activities on soil erosion through land use change can be investigated.
Previous studies on global soil erosion calculated the global R factor based on the total annual precipitation (Renard and Freimund, 1994).This method is different from the method presented in the original RUSLE model (Renard et al., 1997), which is mainly based on 30 min precipitation intensity.The reason for the method of Renard and Freimund (1994) is the lack of high resolution precipitation intensity on a global scale.However, high resolution precipitation intensity is an important explaining parameter of the R factor and therefore the applicability of the method of Renard and Freimund (1994) is limited.
The overall objective of our study is to extend the applicability of the RUSLE model to a coarse resolution at global scale, in order to make the model compatible with ESMs.This would enable future studies on the effects of soil erosion for the past, current and future climate.To this end, we develop generally applicable methods that improve the estimation of slope and climatic factors from coarse resolution global data sets.These methods should not only be applicable across agricultural areas as in the studies of Van Oost et al. (2007) and Doetterl et al. (2012) but also across nonagricultural areas.We adjust the S factor to the coarse resolution of the global scale based on the scaling of slope according to the fractal method.The adjustment of the R factor to the global scale is based on globally applicable regression equations.We derived these regression equations for different climate zones based on parameters for precipitation, elevation and the simple precipitation intensity.This approach is validated using several high resolution data sets on the R factor.Finally, the effects of these adjustments to both factors on global soil erosion rates are investigated separately and tested against independent estimates of soil erosion from high resolution and high precision data sets of Europe and the USA.

Scaling slope according to the fractal method
The topographical factors of RUSLE are the slope steepness factor (S) and a slope length factor (L).The S factor is generally computed by the continuous function of Nearing (1997): (1) And the L factor is computed according to Renard et al. (1997): where in which θ is the slope and l is the slope length in metres.
As seen in Eqs.
(1)-(3), slope is a crucial parameter and thus an accurate estimation is essential in deriving accurate estimates of the L and S factors and soil erosion rates.For an accurate estimation of the slope, input elevation data from digital elevation models (DEMs) should capture the detailed spatial variability in elevation.However, global DEMs are often too coarse to capture the detailed topography because of the surface smoothening effect.To account for this problem it is assumed that topography is fractal.Following Klinkenberg and Goodchild (1992) and Zhang et al. (1999), slope can be expressed as a function of the spatial scale by applying the variogram equation.The variogram equation is used to approximate the fractal dimension of topography and is expressed as follows: so that where Z p and Z q are the elevations at points p and q, d pq is the distance between p and q, k is a constant, α = k 0.5 , and D is the fractal dimension.Because the left side of Eq. ( 5) represents the slope, it can be assumed that the slope (θ ) is related to the spatial scale or the grid size (d) in This result implies that by calculating the fractal properties (D and α) Eq. ( 6) can be used to calculate slope at any specified d.The local fractal dimension (D) describes the roughness of the topography while the local value of α is related to the concept of lacunarity, which is a measure of the size of "gaps" (valleys and plains) in the topography (Zhang et al., 2002).To estimate the spatial variations of D and α, Zhang et al. (1999) proposed to relate these parameters to the standard deviation of elevation.Hereby it is assumed that the standard deviation of elevation does not change much with the DEM resolution.D is then calculated as a function of the standard deviation (σ ) in a 3 pixel × 3 pixel moving window, as proposed by Zhang et al. (1999): To estimate α we used the modified approach by Pradhan et al. (2006).They derived α directly from the steepest slope in a 3 pixel × 3 pixel moving window, called α steepest in the following.Having obtained α steepest and D from a grid at a given resolution, the scaled slope (θ scaled ) for a target grid resolution (d scaled ) is obtained by Pradhan et al. (2006) also showed that in their case study the ideal target resolution for downscaling slope was 150 m.This is due to the breakdown of the unifractal concept at very fine scales, which was shown to happen at a scale of 50 m.Altogether, this fractal method shows that a high resolution slope can be obtained from a low resolution DEM as is needed by the RUSLE model.

Application of the fractal method on global scale
In this study, we investigate the performance of the fractal method on a global scale using different global DEMs as a starting point.The target resolution of downscaling is put to 150 m (about 5 arcsec) according to Pradhan et al. (2006).It should be noted that the spatial scale on which the original RUSLE and USLE models are operating is usually between 10 and 100 m, which indicates that the 150 m target resolution may be still too coarse for a correct representation of slope.The DEMs that are used here are given in Table 1.
As reported in previous studies (Zhang et al., 1999;Chang and Tsai, 1991;Zhang and Montgomery, 1994), the average slope decreases with decreasing DEM resolution.This confirms the expectation of loss of detail in topography at lower DEM resolutions.A large difference is found between the unscaled global average slope from the 5 arcmin and the 30 arcsec DEMs, which is in the order of 0.017 m m −1 or 74 % (Table 2).After applying the fractal method, the scaled slopes at 150 m target resolution from all DEMs increased significantly compared to the unscaled slopes (Fig. 1).However, there is still a difference of about 0.05 m m −1 or 8.5 % between the scaled slopes from the 5 arcmin and the 30 arcsec DEMs (Table 2).This difference can be attributed to several factors.One factor could be the underlying assumption that the standard deviation of elevation (σ ) is independent of the DEM resolution.Although σ does not change much when considering different resolutions, there is still a general decrease in mean global σ when going from the 5 arcmin to the 30 arcsec DEM (Table 2).Due to the dependence of the fractal dimension (D) on σ (Zhang et al., 1999), a decrease of σ leads to a decrease in D and therefore an increase in the scaled slope.Other factors that could play a role here are the dependence of α steepest on the steepest slope, and the breakdown of the fractal method at certain scales and in certain environments.Zhang et al. (1999) mentioned that the scaling properties of slope are affected in very coarse resolution DEMs if σ changes considerably.On the other hand, Pradhan et al. (2006) mentioned the breakdown of the fractal method at very fine scales.This can indicate that the 150 m target resolution is not appropriate for some topographically complex regions in the world or, as addressed by Zhang et al. (1999), the DEMs used in this study are too coarse to scale down the slope to 150 m accurately for these regions.
After applying the fractal method on a 30 arcsec resolution DEM, the scaled slope shows a clear increase in detail, while the unscaled slope shows a strong smoothening effect (Fig. 2a, b).It is found that, after scaling, the slope values range from 0 to 85 • and are less than 2 • in 80 % of the area.In contrast, all slope values are less than 45 • and range between 0 and 2 • in 89 % of this area when slope is computed directly from the 30 arcsec DEM.
The scaled slope from the 30 arcsec DEM will be used in this study to estimate the global soil erosion rates by the RUSLE model.

Adjustment of the rainfall erosivity factor
3.1 The approach by Renard and Freimund (1994) Rainfall erosivity (R factor) is described by Hudson (1971) and Wischmeier and Smith (1978) as the result of the transfer of kinetic energy of raindrops to the soil surface.This causes a detachment of soil and the downslope transport of the soil particles, depending on the amount of energy, rainfall intensity, soil type and cover, topography and management (Da Silva, 2004).The original method of calculating erosivity is described by Wischmeier and Smith (1978) and Renard et al. (1997) as where n is the number of years of records, m j is the number of storms of a given year j , and EI 30 is the rainfall erosivity index of a storm k.The event's rainfall erosivity index EI 30 (MJ mm ha −1 h −1 ) is defined as where e r and v r are, respectively, the unit rainfall energy (MJ ha −1 mm −1 ) and the rainfall depth (mm) during a time period r, and I 30 is the maximum rainfall intensity during a  time period of 30 min (mm h −1 ).The unit rainfall energy, e r , is calculated for each time period as where i r is the rainfall intensity during the time period (mm h −1 ).The information needed to calculate the R factor according to the method of Wischmeier and Smith (1978) is difficult to obtain on a large spatial scale or in remote areas.Therefore, different studies have been done on deriving regression equations for the R factor (Angulo-Martinez et al., 2009;Meusburger et al., 2012;Goovaerts, 1999;Diodato and Bellocchi, 2010).Most of these studies, however, concentrate on a specific area and can therefore not be implemented on the global scale.Studies on global soil erosion estimation by the RUSLE model or a modified version of it (Doetterl et al., 2012;Van Oost et al., 2007;Montgomery, 2007;Yang et al., 2003) have all used the method of Renard and Freimund (1994).Renard and Freimund (1994)  To test how this method performs globally, we calculated the R factor according to the method of Renard and Freimund (1994) (Eq. 12) first.Here we used the 0.25 • resolution annual precipitation data from the Global Precipitation Climatology Centre (GPCC) product (Table 1).Then, we selected three regions to validate the resulting R values and their variability: the USA (EPA, 2001), Switzerland (Meusburger et al., 2012), and the Ebro Basin in Spain (Angulo-Martinez et al., 2009).For these regions, high resolution ero- sivity data are available from pluviographic data of local meteorological stations across the whole region.
Figure 3 shows that the R values computed with the Renard and Freimund (1994) method strongly overestimate R when compared to the high resolution R data of the selected regions.For the USA the R factor of Renard and Freimund (1994) shows an overall overestimation for the western USA and for a large part of the eastern USA when compared to the high resolution R factor (Table 7, Fig. 3a).In particular, a strong overestimation is seen for the northwest coast of the USA.This region is known to have complex rainfall patterns due to the presence of mountains and high local precipitation intensities with frequent snow fall (Cooper, 2011).It should be noted that the USA is not the best suited case study for testing the R values computed with the Renard and Freimund (1994) method, as this method is based on climate data from stations in the USA.The available high resolution or observed data on the R factor from Switzerland and the Ebro Basin are better suited for an independent validation.
For Switzerland, which has a complex precipitation variability influenced by the relief of the Alps (Meusburger et al., 2012), the R factor of Renard and Freimund (1994) shows a strong overall overestimation when compared to the high resolution R values (Table 7, Fig. 3b).For the Ebro Basin, located in Spain, the observed R data were available for the period 1997-2006 from Angulo-Martinez et al. (2009).Also here the method of Renard and Freimund (1994) overestimates the R factor and is not able to reproduce the high spatial variability of the R data (Table 7, Fig. 3c).

The linear multiple regression approach using environmental factors
To better represent the R factor on a global scale, the R estimation was based on the updated Köppen-Geiger climate classification (Table 3, Fig. 4).The Köppen-Geiger climate classification is a global climate classification and is based on the vegetation distribution connected to annual cycles of precipitation and temperature (Lohmann et al., 1993).The reason for this approach is that this classification system includes annual cycles of precipitation and is thus indirectly related to precipitation intensity.Based on this, it is possible to derive regression equations for the R factor that are applicable for each individual climate zone of the classification.This provides a basis to calculate the R factor with coarse resolution data on a global scale.-frost T hot < −0 * MAP: mean annual precipitation, MAT: mean annual temperature, T hot : temperature of the hottest month, T cold : temperature of the coldest month, T mon10 : number of months where the temperature is above 10, P dry : precipitation of the driest month, P sdry : precipitation of the driest month in summer, P wdry : precipitation of the driest month in winter, P swet : precipitation of the wettest month in summer, P wwet : precipitation of the wettest month in winter, P threshold : varies according to the following rules (if 70 % of MAP occurs in winter then P threshold = 2× MAT, if 70 % of MAP occurs in summer then P threshold = 2× MAT + 28, otherwise P threshold = 2× MAT + 14).Summer (winter) is defined as the warmer (cooler) 6-month period of AMJJAS (ONDJFM).
As a basis for deriving the regression equations for the R factor we used high resolution R maps of the USA from the EPA (2001).The USA covers most of the world's climate zones and is also the largest region with available high resolution R data.Linear multiple regression was used to adjust R: where X is the independent explanatory variable, j is the number of explanatory variables, β is a constant and ε is the residual.
The regression operates on one or more of the following parameters (X j ): total annual precipitation (GPCC 0.25 • product), mean elevation (ETOPO 5 DEM), and the simple precipitation intensity index, SDII.It should be mentioned that the SDII was only available on a very coarse resolution of 2.5 • for certain regions on Earth, such as parts of Europe and the USA.The SDII is calculated as the daily precipitation amount on wet days (≥ 1 mm) in a certain time period divided by the number of wet days in that period.Previous studies that performed regression of R showed that precipitation and elevation were in most cases the only explanatory variables (Meusburger et al., 2012;Mikhailova et al., 1997;Goovaerts, 1999;Diodato and Bellocchi, 2010;Angulo-Martinez et al., 2009).Here, we added to the regression the SDII as it is a simple representation of precipitation intensity, which is an important explaining variable of the R factor.The precipitation and SDII data sets were rescaled to a 5 arcmin resolution (corresponding to 0.0833 • ) to match the Köppen-Geiger climate classification data that was available at the resolution of 6 arcmin (corresponding to 0.1 • ).Furthermore, high resolution erosivity data from Switzerland (Meusburger et al., 2012) and annual precipitation from the GPCC 0.5 • product were used to derive the regression equations for the R factor for the polar (E) climate zones.These climate zones are not present in the USA.For the rest of the climate zones that are not present in the USA it was difficult to obtain high resolution erosivity data.Therefore, we maintained the method of Renard and Freimund (1994) for those climate zones to calculate erosivity.Also, we kept the R factor of the Renard and Freimund (1994) method if no clear improvement of the R factor was found when using the new regression equations for a specific climate zone.Here, we mainly used the r 2 combined with the residual standard error to evaluate if the new regression equations  showed a clear improvement in the R factor.The Renard and Freimund (1994) R factors where kept for the hot arid climate zone (BWh) and the temperate climate zone with a hot summer (Csa) in the USA.These are just two climate zones out of the 17 evaluated ones, which show that the Renard and Freimund method performs as good as or slightly better than the regression method.All data sets for deriving the R factor are described in Table 1.

Application of the linear multiple regression method on a global scale
Tables 4 and 5 show the resulting regression equations for climate zones for which we found initially a low correlation between the R values calculated by the method of Renard and Freimund (1994) and the high resolution R values from the EPA (2001) and Meusburger et al. (2012).Figure 5 shows for each addressed climate zone how the method of Renard and Freimund (1994) and the new regression equations compare to the high resolution R of the USA.For the cold climate zones with a dry summer (Ds), the new regression equations show only a slight improvement as compared to the method of Renard and Freimund (1994).Also for the polar climate zones (E) the new regression equations still show a significant bias.However, they perform much better compared to the method of Renard and Freimund (1994).For most of the addressed climate zones the SDII explains a large part of the variability in the R factor.The elevation plays a smaller role here.Elevation can be an important explaining variable in regions with a high elevation variability, which then affects the precipitation intensity.
From Tables 4 and 6 it can be concluded that the R factor in climate zones without a dry season (f) can be easily explained by the total annual precipitation and the SDII.Dry climate zones, especially dry summer climate zones, showed a weaker correlation.This is most likely due to the fact that the SDII is too coarse to explain the variability in the low precipitation intensity in the summer.It is also interesting to see that even though the SDII was derived from a very coarse resolution data set, it turned out to be still important for deriving more accurate R values.
We also show for each addressed climate zone a comparison of the newly computed average R factor with the average high resolution R factor, and the uncertainty range (Table 6).The uncertainty range was computed by taking into account the standard deviation of each of the parameters in the regression equations.As mentioned before, the polar climate zones showed the largest uncertainty range.The new regression equations significantly improved the R values and spatial variability in the western USA and lead to an average R factor that was closer to the data mean (Table 7, Fig. 6a).Although the new regression equations show a bias for the polar climate zones (the minimum and maximum R values are not captured), the resulting mean R values for Switzerland show a strong improvement (Table 7, Fig. 6b).
Furthermore, the variability in the estimated R factor compares well with the variability of the high resolution R factor.It should be noted that Switzerland is not an independent case study for the polar climate zones, as the high resolution R values from this case study were used in our regression analysis.However, the Ebro Basin case study confirms the strong improvement for the polar climate zones (Fig. 6c).As the high resolution R values of the USA and Switzerland were used to derive the regression equations, the third case study, the Ebro Basin in Spain, provided an important independent validation.For the Ebro Basin, the new regression equations not only improve the overall mean but also capture the minimum R values better.This resulted in an improved representation of the R variability (Table 7, Fig. 6c).In Fig. 6c, however, there is a clear pattern separation in the newly computed R values, which is due to the fact that the SDII data are not available for part of the Ebro Basin.As mentioned before, SDII is an important explaining parameter in the regression equations for most of the addressed climate zones.
Figure 7a shows the global patterns of the estimated R factor from the method of Renard and Freimund (1994) and the new regression equations.Figure 7b shows a difference plot between the estimated R factor with the method of Renard and Freimund (1994) and the R factor estimated with the new regression equations.The new regression equations significantly reduced the R values in most regions.However, the tropical regions still show unrealistic high R values (maximum R values go up to 1 × 10 5 MJ mm ha −1 h −1 year −1 ).This is because the R factor was not adjusted for the tropical climate zones due to the lack of high resolution R data.Renard and Freimund (1994) and the new modelled R values (MJ mm ha −1 h −1 year −1 ), where blue colours indicate lower R values by Renard and Freimund (1994) compared to the new modelled R values, while reddish colours indicate higher R values; map resolution is 5 arcmin.Oliveira et al. (2013) found for the R factor in Brazil that the maximum R values for the tropical climate zones reach 22 452 MJ mm ha −1 h −1 year −1 .We find R values in Brazil that exceed this maximum R value found by Oliveira et al. (2013).
Finally, it should be noted that the purpose of the adjusting methods for the S and R factors in this study is to capture more accurately the large-scale mean erosion rates rather than the extremes.Therefore, even though the new regression equations are still not accurate enough for certain climate zones, it is important that the average R factor is represented well.The approach for adjusting the R factor also showed that although there is no high temporal resolution precipitation intensity data available on a global scale, the R factor can still be represented well for most climate zones on a large spatial scale.This can be done by using other parameters, such as elevation, and especially one representative of precipitation intensity, such as the SDII.The SDII played an important role here as it improved the estimation of the R factor significantly, even though data was only available at a very low resolution as compared to the other data sets of precipitation, elevation and climate zone classification.

Computation of the soil erodibility and land cover factors
In the following we demonstrate the consequences of the new parameterizations of the S and R factors for global soil erosion rates.First, we compute the other individual RUSLE factors, soil erodibility (K) and crop cover (C).Estimations of the K factor were based on soil data from the gridded 30 arcsec Global Soil Data set for use in Earth System Mod-  .GSCE is based on the Harmonized World Soil database (HWSD) and various other regional and national soil databases (Shangguan et al., 2014).We used the method of Torri et al. (1997) to estimate the K factor, and gave volcanic soils a K factor of 0.08 t ha h ha −1 MJ −1 mm −1 .This is because these soil types are usually very vulnerable to soil erosion, and the observed K values are beyond the range predicted by the method of Torri et al. (1997) (Van der Knijff et al., 1999).To account for the effect of stoniness on soil erosion we used a combination of the methods by Cerdan et al. ( 2010) and Doetterl et al. (2012), who based their methods on the original method of Poesen et al. (1994).For nonagricultural areas we used the method of Cerdan et al. ( 2010), where they reduced the total erosion by 30 % for areas with a gravel percentage larger or equal to 30 %.For agricultural and grassland areas we used the method of Doetterl et al. (2012), where erosion was reduced by 80 % in areas where the gravel percentage exceeded 12 %.We calculated the C factor according to the method of De Jong et al. (1998), using 0.25 • normalized difference vegetation index (NDVI) and land use data for the year 2002.An important limitation of this method is the fact that in winter the C factor is estimated too high (Van der Knijff et al., 1999).This is because the method does not include the effects of mulch, decaying biomass and other surface cover reducing soil erosion.To prevent the C factor from being too high, maximum C values for forest and grassland of 0.01 and 0.05 for pasture were used.Doetterl et al. (2012) showed that the slope length (L) and support practice (P ) factors do not contribute significantly to the variation in soil erosion at the continental scale to global scale, when compared to the contribution of the other RUSLE factors (S, R and C).However, this does not mean that their influence on erosion should be ignored completely.They may play an important role in local variation of erosion rates.In our erosion calculations we do not include these factors because we have too little or no data of these factors on a global scale.Including them in the calculations would only add an additional large uncertainty to the erosion rates.This would make it more difficult to judge the improvements we made to the S and R factors.

Computation of global soil erosion rates and comparison to empirical databases
We applied the RUSLE model with the settings mentioned in the previous paragraph at a 5 arcmin resolution on a global scale for the present time period (see time resolutions of data sets in Table 1).We calculated global soil erosion rates with four different versions of the RUSLE model: (a) the unadjusted RUSLE, (b) RUSLE with only an adjusted S factor, (c) RUSLE with only an adjusted R factor, and (d) the adjusted RUSLE (all adjustments included).We found a global average soil erosion rate for the adjusted RUSLE of 6.5 t ha −1 year −1 (Fig. 8a).When including the uncertainty arising from applying the linear multiple regression method, the mean global soil erosion rate differs between 5.3 and 15 t ha −1 year −1 .Furthermore, the RUSLE version with only an adjusted S factor shows the highest average global soil erosion rate, while the lowest rate is found for the RUSLE version with only the adjusted R factor (Table 8). Figure 8c shows the difference between the erosion rates of the S-adjusted RUSLE and the unadjusted RUSLE versions.The erosion rates are in general increased here and mostly pronounced in mountainous regions.This feature is "dampened" when adjusting the R factor.The difference between the R-adjusted RUSLE and unadjusted RUSLE versions (Fig. 8d) shows that the erosion rates are overall decreased in regions where the adjustments are made.When the erosion rates of the unadjusted RUSLE model are subtracted from the fully adjusted RUSLE model (Fig. 8b), we find that erosion rates are slightly decreased in areas where the R factor is adjusted.However, for the tropics an increase in erosion rates is found in the fully adjusted RUSLE due to the lack of adjusting the R factor there.This indicates that these two factors balance each other, and that it is important to have a correct representation of all the RUSLE factors on a global scale in order to predict reliable erosion rates.
In this study the K and C factors are not tested and adjusted for a coarse resolution at global scale and thus validation with existing empirical databases on soil erosion is not fully justified.However, to test if the global erosion rates are in an acceptable range, they are compared to erosion estimates from the NRI database for the USA and erosion estimates from the study of Cerdan et al. (2010) for Europe.These are to our knowledge the only large-scale high resolution empirical databases on soil erosion.
The NRI database contains USLE erosion estimates for the year 1997, which are available at the Hydrologic Unit Code 4 (HUC4) watershed level.We aggregated the resulting erosion rates from the adjusted and unadjusted RUSLE models to the HUC4 watershed level.The results show that the average erosion rates from the adjusted RUSLE model come closer to that of the NRI database (Table 9, Fig. 9a).However, the maximum average HUC4 soil erosion rate from the adjusted RUSLE is somewhat higher compared to the NRI database.From these results we can conclude that the erosion rates of the adjusted RUSLE fall in the range of observed values but that there are still some local overestimations.Some of these overestimations can be found in the south-west of the USA where the adjusted RUSLE shows a slightly worse performance compared to the unadjusted RUSLE.The R factor in this region was not changed as it was already estimated well by the method of Renard and Freimund (1994), however, the S factor increased due to the hilly terrain.Without adjusting the other RUSLE factors (K and C), this resulted in an overall increase in soil erosion rates.This indicates that the other RUSLE factors may play an important role in this region.Furthermore, we see that along the west coast of the USA the erosion values are not much improved with the adjusted RUSLE model.This is mainly because some climate zones such as the temperate climate zone with a dry and warm summer (Csb) prevail in this region, for which the R factor is still difficult to estimate in a correct way (Table 4).
For Europe, Cerdan et al. ( 2010) used an extensive database of measured erosion rates on plots under natural rainfall.They extrapolated measured erosion rates to all of Europe (European Union area) and adjusted them with a topographic correction.This correction was based on the L and S factors of the RUSLE model.They also applied a correction to account for soil stoniness.For comparison, the soil erosion rates from Cerdan et al. ( 2010) and the RUSLE estimates in our study are aggregated at country level.The performance of the adjusted RUSLE model was not as good for Europe as compared to the USA.This is not surprising as the RUSLE model is based on soil erosion data of the USA.However, also on the European scale the adjusted RUSLE model performed better than the unadjusted RUSLE model (Table 9, Fig. 9b).In particular, the large erosion rates in the south of Europe as observed in the results of the unadjusted RUSLE model are less extreme in the adjusted RUSLE model.Still, the overall average erosion rate for Europe is overestimated by approximately 2 times (Table 9).
The biases in erosion rates as seen for the south-west of the USA and southern Europe can be attributed to several factors.As mentioned before, the other RUSLE factors (K and C) and the way they interact with the R and S factors are not adjusted to the coarse resolution at global scale.We found no clear signal for the land cover types with which the adjusted RUSLE performs better or worse.In general, we can see that the adjusted RUSLE model still overestimates erosion rates for most land cover types.A short analysis for Europe showed that the largest biases are found for shrubs and the lowest for grassland.However, a more explicit analysis is needed to find out how we can improve the contribution of land cover and land use to erosion rates in the RUSLE model.Explicitly including the interaction between the C and R factor on a monthly timescale could be crucial.This is very important for example in areas with agriculture and areas with a strong seasonal character.Another aspect related to improving the C factor is looking at the location of land use in a certain grid cell.If the land use in a grid cell is located on steep slopes, the resulting erosion in that grid cell would be higher than when it was located in the flatter areas.In this study, however, only mean fractions of land cover and the NDVI are used for each grid cell.This can lead to possible biases in the resulting erosion rates.
Furthermore, land management is not accounted for in this study, which could introduce an important systematic bias in the soil erosion rates especially for agricultural areas.Land management is represented by the P factor in the original USLE; however, it is partly also incorporated in the C factor for agricultural land use through plant residues, cover crops and tillage.A limitation of the NDVI approach to estimate the C factor lies therefore in the inability to estimate this land management effect.Applying this method also limits the interaction between the R and C factors on a monthly to seasonal scale, because this interaction is partly based on land management.
Furthermore, uncertainties in the coarse resolution land cover/land use, soil and precipitation data sets that are not accounted for can lead to the model biases.Also, better adjustment of the R factor for climate zones such as the polar climates could help improve the overall results.Some biases in the erosion rates can also be attributed to the fact that stepped relief, where flat plateaus are separated by steep slopes, is not well captured by the 150 m target resolution used in the fractal method to scale slope.In this way erosion would be overestimated in these areas.Finally, errors and limitations in the observational data sets can also contribute to the dif-ferences between model and observations.The study of Cerdan et al. ( 2010) on Europe, for example, used extrapolation of local erosion data to larger areas, which could introduce some biases.Also, the underlying studies on measured erosion rates used different erosion measuring techniques that can be linked to different observational errors.

Conclusions
In this study we introduced specific methods to adjust the topographical and rainfall erosivity factors to improve the application of the RUSLE model on global scale, using coarse resolution input data.
Our results show that the fractal method by Zhang et al. (1999) and Pradhan et al. (2006) can be applied on coarse resolution DEMs to improve the resulting slope.Although the slope representation improved after applying this method, the results still show a slight dependence on the original grid resolution.This is attributable to several factors such as the underlying assumption that the standard deviation of eleva-tion (σ ) is independent of the DEM resolution and to the breakdown of the fractal method at certain scales.
We compared the rainfall erosivity calculated by the method of Renard and Freimund (1994) to available high resolution or observed erosivity data of the USA, Switzerland and the Ebro Basin.We find that this method results in overall significant biases in erosivity.Therefore, we implemented a linear multiple regression method to adjust erosivity for climate zones of the Köppen-Geiger climate classification system in the USA.Using precipitation, elevation and the simple precipitation intensity index as explaining parameters, the resulting adjusted erosivity compares much better to the observed erosivity data for the USA, Switzerland and the Ebro Basin.Not only are the mean values improved but also the spatial variability in erosivity.It was surprising to notice that using the rather coarse resolution simple precipitation intensity index in the regression analysis made it possible to explain much of the variability in erosivity.This, once more, underpins the importance of precipitation intensity in erosivity estimation.
After calculating the newly adjusted erosivity on a global scale, it is apparent that the tropical climate zones, for which erosivity was not adjusted, show strong overestimations in some areas.This shows that adjusting erosivity for the tropical climate zones should be the next step.The challenge is to find enough reliable long-term and high resolution erosivity data for those regions.
To investigate how the adjusted topographical and rainfall erosivity factors affect the global soil erosion rates, we applied the adjusted RUSLE model on a global scale.We found an average global soil erosion rate of 6.5 t ha −1 year −1 .It is, however, difficult to provide accurate uncertainty estimates to these global erosion rates and to provide a good validation with observations.This is due to lack of high resolution data on other individual RUSLE factors such as the land cover, soil erodibility, slope length and support practice.These RUSLE factors are therefore not adjusted for application at coarse resolution on a global scale.We argue that it is important to focus on adjusting the other RUSLE factors for an improved application of the RUSLE model on global scale.The next step would be to better capture the anthropogenic contribution to global soil erosion.This can be done by adjusting first of all the land cover factor to a coarse resolution application and focusing on the interaction of this factor with rainfall erosivity on a monthly to seasonal basis.This is important because the land cover factor has strong interactions with the rainfall erosivity factor and includes the effect of human activities on erosion through agricultural activities and land management.
To test if the soil erosion rates from the adjusted RUSLE model are in a realistic range, we compared the results to the USLE erosion estimates for the USA from the NRI database and the erosion estimates for Europe from the study of Cerdan et al. (2010).The adjusted RUSLE soil erosion rates, which we aggregated to the watershed level, show a better comparison with the NRI USLE estimates than the unadjusted RUSLE erosion rates.For Europe, the comparison of the adjusted RUSLE soil erosion rates to the study of Cerdan et al. (2010) were not as good as for the USA.This is not surprising due to the fact that the parameterizations of the RUSLE model are based on soil erosion data of the USA.However, also for Europe, the adjusted RUSLE model performs better than the unadjusted RUSLE model.
We find overestimations by the adjusted RUSLE model for hilly regions along the west coast of the USA and for southern of Europe.We argue that, besides the reasons mentioned before, these biases are due to the fact that the topographical detail may not be enough in some regions to capture the true variability in soil erosion effects by topography.Also, erosivity could not be adjusted for some climate zones that are not present in the USA or Switzerland and needs to be further improved for climate zones such as the polar climate zones.
We conclude that even though there is still much improvement possible in the RUSLE model with respect to topography and erosivity, the methods proposed in this study seem to be promising tools for improving the global applicability of the model.A globally applicable version of the RUSLE model, together with data on environmental factors from ESMs, can be a basis for future studies on accurate soil erosion rates for past, current and future scenarios.

V.
Naipal et al.: Improving the global applicability of the RUSLE model

Figure 1 .
Figure 1.Global average unscaled slope estimated from different coarse resolution digital elevation models (DEMs) as function of their resolution (blue), and global average scaled slope from the same DEMs as function of their resolution (red).

Figure 2 .
Figure 2. (a) A global map of the scaled slope derived from the 30 arcsec DEM using a target resolution of 150 m.(b) A global map showing the difference between the unscaled and scaled slopes (in degrees), where blue colours show an underestimation by the unscaled slope when compared to the scaled slope and reddish colours show and overestimation.
related the R factor to the total annual precipitation based on erosivity data available for 155 stations in the USA, shown in the following equations: R = 0.0483 × P 1.61 , P ≤ 850 mm R = 587.8− 1.219 × P + 0.004105 × P 2 , P >850 mm.(12)

Figure 3 .
Figure 3. Spatial difference plots showing the difference between the high resolution R values and R values calculated with the method of Renard and Freimund (1994) for (a) the USA, (b) Switzerland and (c) the Ebro Basin in Spain; in panels (a) and (b) the blue colours show an underestimation of the calculated R factor when compared to the high resolution R values, while the red colours show an overestimation; the Ebro Basin serves here as an independent validation set and it has two graphs: (c1) a spatial plot of erosivity according to Renard and Freimund (1994) and (c2) the high resolution R values from Angulo-Martinez et al. (2009) (all values in the graphs are in MJ mm ha −1 h −1 year −1 ).

Table 4 .
Linear multiple regression equations for different climate zones, relating high resolution R factor from the USA with one or more significant parameters: annual total mean precipitation, P (mm), mean elevation, z (m), and the simple precipitation intensity index, SDII (mm day −1

Figure 5 .
Figure 5.Comparison of high resolution R factor data and predicted R values from (1) theRenard and Freimund (1994) method and(2) the new regression equations, for various climate zones; the red line is the 1-to-1 line and does not appear in some graphs because predicted R values are overestimated.

Figure 6 .
Figure 6.Spatial difference plots showing the difference between the high resolution R values and R values calculated with the new regression equations for (a) the USA, (b) Switzerland and (c) the Ebro Basin in Spain; in panels (a) and (b) the blue colours show an underestimation of the calculated R values when compared to the high resolution R values, while the red colours show an overestimation; the Ebro Basin serves here as an independent validation set and it has two graphs: (c1) a spatial plot of the R factor according to the new regression equations, and (c2) the high resolution R values from Angulo-Martinez et al. (2009) (all values in the graphs are in MJ mm ha −1 h −1 year −1 ).

Figure 7 .
Figure 7. (a) Global distribution of the new modelled R values according to the new regression equations; and (b) a difference map between R values calculated according to the method ofRenard and Freimund (1994) and the new modelled R values (MJ mm ha −1 h −1 year −1 ), where blue colours indicate lower R values byRenard and Freimund (1994) compared to the new modelled R values, while reddish colours indicate higher R values; map resolution is 5 arcmin.

Figure 8 .
Figure 8.(a) Global yearly averaged erosion rates according to the fully adjusted RUSLE model; (b) a difference map between the fully adjusted and unadjusted RUSLE model; (c) a difference map between the adjusted S-RUSLE model and the unadjusted RUSLE model; (d) a difference map between the adjusted R-RUSLE model and the unadjusted RUSLE model.In panels (b), (c) and (d) the reddish colours show an overestimation of the adjusted RUSLE model and yellow to bluish colours show an underestimation (resolution of all maps is 5 arcmin and all units are in t ha −1 year −1 ).

Figure 9 .
Figure 9. (Top) Difference plots between soil erosion estimates from the NRI database for the USA and estimates of (a) the unadjusted RUSLE model, and of (b) the adjusted RUSLE model, all aggregated at HUC4 watershed level.(Bottom) Difference plots between soil erosion estimates from the database of Cerdan et al. (2010) for Europe and estimates of (c) the unadjusted RUSLE model and of (d) the adjusted RUSLE model, all aggregated at country level.Reddish colours represent an overestimation (t ha −1 year −1 ) while the bluish colours represent and underestimation (t ha −1 year −1 ) compared to the erosion values from the databases.

Table 1 .
List of data sets used in this study.

Table 2 .
Fractal parameters and the resulting mean global slopes before and after applying the fractal method on the different DEMs.Increase of slope means the increase of the average global slope of a DEM after applying the fractal method; difference after scaling = www.geosci-model-dev.net/8/2893/2015/Geosci.Model Dev., 8, 2893-2913, 2015

Table 5 .
Linear multiple regression equations for different climate zones for regions that have no data on the simple precipitation intensity index, SDII (mm day −1 ).The regression equations relate high resolution erosivity from the USA to the annual total mean precipitation, P (mm), and/or the mean elevation, z (m).

Table 6 .
Mean high resolution R values (MJ mm ha −1 h −1 year −1 ) from the USA and Switzerland and mean modelled R values with uncertainty range for each addressed climate zone.

Table 7 .
Statistics of the comparison of high resolution R

Table 8 .
Comparison of the global erosion rates (t ha −1 year −1 ) and percentiles between different versions of the RUSLE model.

Table 9 .
Statistics of the observed and modelled erosion rates from the unadjusted and adjusted versions of the RUSLE for the USA and Europe (t ha −1 year −1 ).