# Global rainfall erosivity assessment based on high-temporal resolution rainfall records

Panos Panagos, Pasquale Borrelli, Katrin Meusburger, Bofu Yu, Andreas Klik, Kyoung Jae Lim, Jae Yang, Jinren Ni, Chiyuan Miao, Nabansu Chattopadhyay, et al.

► To cite this version:

Panos Panagos, Pasquale Borrelli, Katrin Meusburger, Bofu Yu, Andreas Klik, et al. Global rainfall erosivity assessment based on high-temporal resolution rainfall records. Scientific Reports, 2017, 7, pp.4175. 10.1038/s41598-017-04282-8. hal-03226776

HAL Id: hal-03226776

https://hal.science/hal-03226776

Submitted on 16 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

# Global rainfall erosivity assessment based on high-temporal resolution rainfall records

Received: 9 February 2017

Accepted: 2 June 2017

Published online: 23 June 2017

Panos Panagos¹, Pasquale Borrelli², Katrin Meusburger², Bofu Yu³, Andreas Klik⁴, Kyoung Jae Lim⁵, Jae E. Yang⁵, Jinren Ni⁶, Chiyuan Miao⁷, Nabansu Chattopadhyay⁸, Seyed Hamidreza Sadeghi⁹, Zeinab Hazbavi⁹, Mohsen Zabihi⁹, Gennady A. Larionov¹⁰, Sergey F. Krasnov¹⁰, Andrey V. Gorobets¹⁰, Yoav Levi¹¹, Gunay Erpul¹², Christian Birkel¹³, Natalia Hoyos¹⁴, Victoria Naipal¹⁵, Paulo Tarso S. Oliveira¹⁶, Carlos A. Bonilla¹⁷, Mohamed Meddi¹⁸, Werner Nel¹⁹, Hassan Al Dashti²⁰, Martino Boni¹, Nazzareno Diodato²¹, Kristof Van Oost²², Mark Nearing²³ & Cristiano Ballabio¹

The exposure of the Earth's surface to the energetic input of rainfall is one of the key factors controlling water erosion. While water erosion is identified as the most serious cause of soil degradation globally, global patterns of rainfall erosivity remain poorly quantified and estimates have large uncertainties. This hampers the implementation of effective soil degradation mitigation and restoration strategies. Quantifying rainfall erosivity is challenging as it requires high temporal resolution (<30 min) and high fidelity rainfall recordings. We present the results of an extensive global data collection effort whereby we estimated rainfall erosivity for 3,625 stations covering 63 countries. This first ever Global Rainfall Erosivity Database was used to develop a global erosivity map at 30 arc-seconds (~1 km) based on a Gaussian Process Regression (GPR). Globally, the mean rainfall erosivity was estimated to be 2,190 MJ mm ha⁻¹ h⁻¹ yr⁻¹, with the highest values in South America and the Caribbean countries, Central east Africa and South east Asia. The lowest values are mainly found in Canada, the Russian Federation, Northern Europe, Northern Africa and the Middle East. The tropical climate zone has the highest mean rainfall erosivity followed by the temperate whereas the lowest mean was estimated in the cold climate zone.

Given the growing concerns about climate change, climatic data is particularly important for the scientific community and society in general, as decisions of individuals, business and governments are dependent on available meteorological data¹. At present, a large number of the large-scale precipitation datasets are publicly available, although formats and completeness of the records vary widely². Heavy rainfall and extreme events are of major importance for climate change, economy and society³. However, extreme events are typically rare events of short

¹European Commission, Joint Research Centre, I-21027, Ispra (VA), Italy. ²Environmental Geosciences, University of Basel, Basel, Switzerland. ³School of Engineering, Griffith University, Nathan, Australia. ⁴BOKU, University of Natural Resources and Life Sciences, Vienna, Austria. ⁵Kangwon National University, Chuncheon-si, Gangwon-do, South Korea. ⁶College of Environmental Sciences and Engineering, Peking University, Beijing, P.R. China. ⁷College of Global Change and Earth System Science, Beijing Normal University, Beijing, P.R. China. ⁸India Meteorological Department, Pune, India. ⁹Faculty of Natural Resources, Tarbiat Modares University, Jalal, Iran. ¹⁰Faculty of Geography, Lomonosov Moscow State University, Moscow, Russian Federation. ¹¹Israel Meteorological Service, Beit Dagan, Israel. ¹²Faculty of Agriculture - Soil Science Department, Ankara University, Ankara, Turkey. ¹³University of Costa Rica, San Jose, Costa Rica. ¹⁴Universidad del Norte, Barranquilla, Colombia. ¹⁵Laboratoire des Sciences du Climat et de l'Environnement, IPSL-LSCE, Gif sur Yvette, France. ¹⁶Federal University of Mato Grosso do Sul, Campo Grande, Brazil. ¹⁷Departamento de Ingeniería Hidráulica y Ambiental, Pontificia Universidad Católica de Chile, Región Metropolitana, Chile. ¹⁸Ecole Nationale Supérieure d'Hydraulique de Blida, Soumaâ, Algeria. ¹⁹Department of Geography and Environmental Science, University of Fort Hare, Alice, South Africa. ²⁰Department of Meteorology, Kuwait, Kuwait. ²¹Met European Research Observatory, Benevento, Italy. ²²Université Catholique de Louvain, Louvain-la-Neuve, Belgium. ²³USDA-ARS, Southwest Watershed Research Center, Tucson, AZ, USA. Correspondence and requests for materials should be addressed to P.P. (email: panos.panagos@ec.europa.eu)
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Global distribution of rainfall erosivity stations (red dots) compiled in the Global Rainfall Erosivity Database (GloREDa); (b) Distribution of rainfall erosivity stations by continent. Maps generated with ESRI ArcGIS ver. 10.4 (http://www.esri.com).</figcaption></figure>

duration and in many regions of the world only limited observational data of sufficient temporal resolution is available to capture them⁴. High temporal resolution rainfall measurements are important in this context, but also have been recognised to be very useful for urban drainage models⁵, climate change modelling⁶, cropping patterns and crop production⁷.

Further, the patterns of heavy and violent rainfall, as captured by the rainfall erosivity factor, influence hydrological and erosive processes and as such are essential for the definition of soil and water conservation practices in the adaptation of agriculture to climate change⁸. Rainfall erosivity is one of the most important input parameters for describing erosive processes and proposing conservation measures by using soil erosion prediction models. Since soil erosion is difficult to measure at large scales, models are required for estimating soil loss by water erosion at regional, national and continental scales. Large scale and global model predictions are of utmost importance, since soil erosion is, in addition to soil sealing, the major threat to soil sustainability and consequently to water- and food security.

As a consequence, recent developments in the soil modelling and climate change communities aim at addressing major scientific gaps in describing key soil processes such as water erosion⁹, based on updated global datasets. Nonetheless, policy initiatives such as the United Nations Convention to Combat Desertification (UNCCD) for zero net land degradation¹⁰, the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES), the Land Degradation and Restoration Assessment¹¹, and the Food Agricultural Organisation (FAO) World Soil Resources¹² still remarked the lack of an updated scientific dataset on global soil erosion.

A vital component of such a global soil erosion map is a spatial assessment of rainfall erosivity. However, methods that estimate erosivity based on annual rainfall data are problematic and highly biased, since rainfall intensity is typically not considered¹³. In order to include rainfall intensity on the calculation of rainfall erosivity, it is necessary to have high temporal resolution rainfall data for long time periods. In this study, we aim to tackle the challenging task of compiling and processing the first global erosivity dataset from long-term high-resolution rainfall data using sub-hourly and hourly pluviographic records. We have used this global erosivity station dataset combined with a set of exhaustive secondary environment variables to generate a global rainfall erosivity map in order to improve our understanding of the global patterns of high intensity rainfall events.

## Results

**Global Rainfall Erosivity Database - GloREDa.** At global scale, this is the first time ever that an erosivity database of such dimension is compiled. The Global Rainfall Erosivity Database, named hereafter as GloREDa, contains erosivity values estimated as R-factors (refer to the method section) from 3,625 stations distributed in 63 countries worldwide. This is the result of an extensive data collection of high temporal resolution rainfall data from the maximum possible number of countries in order to have a representative sample across different climatic and geographic gradients. GloREDa has three components, which are described in the methods section: a) the Rainfall Erosivity database at European Scale (REDES)¹⁴ b) 1,865 stations from 23 countries outside Europe and c) 85 stations collected from a literature review.

The number of GloREDa stations varied greatly among continents (Fig. 1). Europe had the largest contribution to the dataset, with 1,725 stations (48% of total), while South America had the lowest number of stations (141 stations or ~4% of total). Africa has very low density of GloREDa stations (5% of the total). In North America and the Caribbean, we collected erosivity values from 146 stations located in 6 countries (Unites States, Canada, Mexico, Cuba, Jamaica and Costa Rica). Finally, Asia and the Middle East were well represented in GloREDa, with 1,220 stations (34% of the total) distributed in 10 countries including the Siberian part of the Russian Federation (Fig. 1b). The geographic distribution within each continent also differed substantially. For instance, stations in Europe, Oceania and North America covered most of the territory, while those in Africa and South America were largely clustered. However, the stations are well distributed among different erosivity classes (Fig. 2b).
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Global Rainfall Erosivity map (spatial resolution 30 arc-seconds). Erosivity classes correspond to quantiles. Map generated with ESRI ArcGIS ver. 10.4 (http://www.esri.com); (b) number and cumulative percentage of GloREDa stations grouped by erosivity; (c) mean erosivity by continent; (d) mean erosivity by climate zone.</figcaption></figure>

**Global erosivity map.** The Gaussian Process Regression (GPR) model used to interpolate the erosivity (R-factor) point values to a map showed a good performance for the cross-validation dataset [$R^2$ = 0.722, RMSE (Root Mean Square Error) = 1,629 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$]. The annual global erosivity map (Fig. 2) is presented at 30 arc-seconds (~1 km) spatial resolution and subdivided in 10 erosivity classes corresponding to the quantiles. The mean of the global R-factor map is 2,190 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$ with high variability as expressed by the standard deviation of 2,974 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$. The median (50th percentile) of the global erosivity map is 1,150 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$ while 20% of the erosivity values (20th percentile) are lower than 200 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$ and the highest 20% (80th percentile) are higher than 5,200 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$ (Fig. 2). According to the global erosivity map, the highest values are located in south-eastern Asia (Cambodia, Indonesia, Malaysia, the Philippines and Bangladesh), Central Africa (Congo and Cameroon), South America (Brazil, Colombia and Peru), Central America and the Caribbean. The lowest erosivity is mainly located in Siberia, the Middle East, Northern Africa, Canada and Northern Europe. The Polar Regions have been masked out in the global erosivity map.

We found that the spatial patterns of the highest erosivity values (Fig. 2) are coincident with the corresponding patterns of extreme rainfall events reported by Zipser et al.¹⁴. Zipser et al.¹⁴ defined intense storms based on the convective vertical velocity of rain. The authors compiled a 7-year period (1998–2004) database of intense storms and they generated global maps of extreme rainstorm events based on lightning flash, brightness, temperature and noise. According to their study the highest frequency of extreme rainfall events (similar to high annual erosivity values) occurs in the central part of Latin America, Gulf of Mexico, central and western Africa, Madagascar, south-eastern Asia (mainly Bangladesh, south China), Indonesia and North Australia.

**Continental assessments.** At the continental level, South America experiences the highest mean R-factor with 5,874 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$, followed by Africa (3,053 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$), Asia and the Middle East (1,487 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$). In Oceania, the mean R-factor was estimated at 1,675 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$ (Fig. 2c).

Africa exhibits the highest erosivity estimates at the country level; Mauritius and Comoros have the highest worldwide mean annual erosivity values with an R-factor close to 20,000 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$. In Western Africa (Liberia, Sierra Leone and Equatorial Guinea), Central Africa (D.R of Congo, Republic of Congo and Cameroon) and Madagascar mean annual R-factor is higher than 7,000 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$. These patterns agree with those from other continental-scale assessments¹⁵,¹⁶ which indicated highest erosivity values (>10,000 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$) along the Guinea coast of west and central Africa, the Congo basin and Madagascar. Ethiopia and South Africa have mean R-factor values close to 2,500 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$, but the spatial patterns are highly variable with the Ethiopian highlands having extremely high erosivity (>7,000 MJ mm ha$^{-1}$ h$^{-1}$ yr$^{-1}$) while the lowlands
have 3–4 times smaller values. The lowest mean R-factor, with values less than 115 MJ mm ha⁻¹ h⁻¹ yr⁻¹, was estimated for Western Sahara, Libya and Egypt.

Within Asia, the Middle East has the lowest erosivity values, with a mean annual R-factor less than 220 MJ mm ha⁻¹ h⁻¹ yr⁻¹ in Jordan, Saudi Arabia, Kuwait, Syria, Iran and Iraq (Fig. 2). China has a mean value of 1,600 MJ mm ha⁻¹ h⁻¹ yr⁻¹, but exhibits high variability with zero erosivity in the arid north-west areas (Taklimakan desert), and extreme erosivity (>15,000) in the south-eastern coastal zones. Regional studies conducted by Zhu and Yu¹⁷ and Qin et al.¹⁸ show very similar spatial patterns compared to our rainfall erosivity distribution in China. In Japan, the mean annual erosivity was estimated as 4,815 MJ mm ha⁻¹ h⁻¹ yr⁻¹, a value close to the 5,130 MJ mm ha⁻¹ h⁻¹ yr⁻¹ modelled by Shiono et al.¹⁹.

As expected the Siberian part of the Russian Federation and the former Union of Soviet Socialist Republics (Kazakhstan, Turkmenistan and Uzbekistan) have very low mean erosivity values (<250 MJ mm ha⁻¹ h⁻¹ yr⁻¹) due their continental climate. On the contrary, Southeast Asia falls almost completely within the highest erosivity class (>7,400 MJ mm ha⁻¹ h⁻¹ yr⁻¹), in agreement with national assessments for Peninsular Malaysia²⁰. Their erosivity values, generated from pluviographic data range from 7,500 to 20,000 MJ mm ha⁻¹ h⁻¹ yr⁻¹.

In South America, Chile has the lowest R-factor with a mean annual value 1,320 MJ mm ha⁻¹ h⁻¹ yr⁻¹, followed by Argentina (2,232 MJ mm ha⁻¹ h⁻¹ yr⁻¹). The rest of the South American countries have high mean erosivity values (>3,700 MJ mm ha⁻¹ h⁻¹ yr⁻¹), with the highest ones in Brazil, Colombia and Ecuador (>7,000 MJ mm ha⁻¹ h⁻¹ yr⁻¹). The erosivity gradient created by the Andes is clearly visible in the erosivity map. There were few national assessments on rainfall erosivity in south America²¹⁻²⁴. Most of the data used for these studies have been used as input for GloREDa and their spatial patterns are in broad agreement to ours.

In North America and the Caribbean, the mean R-factor is 1,409 MJ mm ha⁻¹ h⁻¹ yr⁻¹ with very low values in Canada and the Northern part of the United States, and extremely high values (>8,000 MJ mm ha⁻¹ h⁻¹ yr⁻¹) along the Gulf of Mexico and the Caribbean countries. The erosivity map for the United States²⁵ also shows high values along the Gulf of Mexico and southern Florida (>8,500 MJ mm ha⁻¹ h⁻¹ yr⁻¹), while overall low values are observed in the Midwestern region (<690 MJ mm ha⁻¹ h⁻¹ yr⁻¹).

In Australia, the mean R-factor is 1,535 MJ mm ha⁻¹ h⁻¹ yr⁻¹ close to 1,767 MJ mm ha⁻¹ h⁻¹ yr⁻¹ estimated by Teng et al.²⁶ based on 11 years (2002–2012) satellite derived Tropical Rainfall Measuring Mission data. In terms of spatial patterns, Teng et al.²⁶ also found maximum erosivity values along the northern and eastern coastal areas (>8,000 MJ mm ha⁻¹ h⁻¹ yr⁻¹), which decreased towards the south-central region (<300 MJ mm ha⁻¹ h⁻¹ yr⁻¹). In New Zealand, the high erosivity values (>4,000 MJ mm ha⁻¹ h⁻¹ yr⁻¹) occur on the west coast of the South Island and decrease towards the east similar to the patterns observed by Klik et al.²⁷ based on 35 weather stations.

The mean erosivity value for Europe was 488 MJ mm ha⁻¹ h⁻¹ yr⁻¹, which is much lower than the one estimated by Panagos et al.²⁸ for the European Union (722 MJ mm ha⁻¹ h⁻¹ yr⁻¹). This is due to the inclusion of European Russia, Ukraine (422 MJ mm ha⁻¹ h⁻¹ yr⁻¹) and Belarus (365 MJ mm ha⁻¹ h⁻¹ yr⁻¹), all of which have low values compared to the other European countries.

**Analysis by Climate zones.** The global rainfall erosivity map was further analysed per climate zones. The updated world Kopper-Geiger climate classification²⁹ is the most widely used and accepted climate map in the scientific community. As expected, tropical climate group showed the highest mean erosivity with 7,104 MJ mm ha⁻¹ h⁻¹ yr⁻¹. Within this group the tropical rainforest (Af) and monsoon (Am) climatic types had the highest mean erosivity and the lowest variability (Fig. 3). Second highest mean erosivity (3,729.3 MJ mm ha⁻¹ h⁻¹ yr⁻¹) occurs in the temperate climate group.

The humid temperate, and temperate with dry winter climate type (Cfa, Cwa), mainly present in the south-eastern United States, eastern Australia and southeast China, have mean erosivity values higher than 4,600 MJ mm ha⁻¹ h⁻¹ yr⁻¹. The Mediterranean (Csa, Csb) and the Oceanic (Cfb) climate zones have mean erosivity values lower than 2,000 MJ mm ha⁻¹ h⁻¹ yr⁻¹ (Fig. 3).

The arid climate group has a relatively low mean erosivity (842 MJ mm ha⁻¹ h⁻¹ yr⁻¹) characterised by the highest spatial variability (e.g. the Cold desert (BWk) type). In this group, the hot desert (BWh) has the largest spatial share (13.9% of global area) with low mean erosivity values (537 MJ mm ha⁻¹ h⁻¹ yr⁻¹). The cold desert climate (BWk), characteristic of northwest China and large areas of Kazakhstan, Turkmenistan, Uzbekistan, North Chile and Argentina, has a very low mean erosivity of 362 MJ mm ha⁻¹ h⁻¹ yr⁻¹. The hot steppe climate (BSh), which is a transition from hot desert to the tropical group (mainly in Africa and India), had medium mean erosivity of 2,371 MJ mm ha⁻¹ h⁻¹ yr⁻¹.

The cold climate group had the lowest mean erosivity, with 493 MJ mm ha⁻¹ h⁻¹ yr⁻¹ whereof the subarctic or boreal climate type (Dsc, Dwc, Dfc), covering major areas of Scandinavia, Siberia and Canada, had minimum mean erosivity values (<285 MJ mm ha⁻¹ h⁻¹ yr⁻¹). By comparison, the climate zones immediately north of hot continental summers (Dfb, Dwb) that cover most of central and eastern Europe, European Russia and the northern United States, have much higher mean erosivity values (526 MJ mm ha⁻¹ h⁻¹ yr⁻¹). The polar areas, mainly located in the Alps, Pyrenees and part of the Tibetan plateau, have a mean erosivity of approximately 990 MJ mm ha⁻¹ h⁻¹ yr⁻¹.

The greatest uncertainty of the global erosivity map is likely related to transition areas between different climatic zones. The different climatic conditions, which result in high variability of rainfall amount, duration, magnitude and intensity, is the main reason for different spatial patterns of erosivity between climatic zones. The standard deviation shows the variability inside the climatic zone (Fig. 3). Moreover, the seasonal variation of climatic conditions play an important role in rainfall erosivity variability.

## Discussion

The increasing availability of rainfall data with high temporal resolution, the growing computing power, and the development of sophisticated geostatistical models, enabled the development of a global rainfall erosivity dataset
<figure><img src="image_4.png" /></figure>

**Figure 3.** R-factor descriptive statistics per Kopper-Geiger climate type. Colour bars are the mean values per climate zone. Error bars represent the standard deviation. Percentages below each main climate category represent its proportion within the study area. Climate zones: Af (tropical rainforest), Am (tropical monsoon), Aw (tropical savannah), BWh (hot desert), BWk (cold desert), BSh (hot steppe), BSk (cold steppe), Csa (dry hot summer), Csb (dry warm summer), Cwa (subtropical dry winter), Cwb (dry winter and dry summer), Cfa (temperate without dry season and hot summer), Cfb (temperate without dry season and warm summer), Cfc (temperate without dry season and cold summer), DSa (cold and dry hot summer), Dsb (cold and dry warm summer), Dsc (cold and dry cold summer), Dwa (cold and dry winter, and hot summer), Dwb (cold and dry winter, and warm summer), Dwc (cold and dry winter, and cold summer), Dwd (cold and dry winter, and very cold winter), Dfa (cold without dry season and hot summer), Dfb (cold without dry season and warm summer), Dfc (cold without dry season and cold summer), Dfd (cold without dry season and very cold winter), E (polar).

at 30 arc-seconds (~1 km) spatial resolution. We acknowledge that this achievement was only feasible through the scientific cooperation between scholars from all over the globe. The global erosivity map was possible thanks to the contribution of data providers (see the long list of meteorological services, organisations, and institutions in the Acknowledgements section), tested methodologies and geostatistical models suitable for such a scale.

**Comparison with past studies.** Compared to previous works on global R-factor estimation³⁰⁻³², our study presents a data-driven approach including measured hourly and sub-hourly data on rainfall intensity for erosivity assessment. In contrast previous global erosivity assessments by Nachtergaele et al.³⁰ and Yang et al.³² that used annual rainfall to estimate erosivity. However, as demonstrated by past literature¹³, ²⁸, ³³ and by the erosivity density statistics, the relation between rainfall amount and erosivity is not straightforward and as shown by Naipal et al.³¹ lead to an overestimation of rainfall erosivity. As a consequence, Naipal et al.³¹, included a simple precipitation intensity index calibrated to high resolution R-factor data of the USA.

In order to compare the three past erosivity maps³⁰⁻³², and the current one, four Ordinary Least Square (OLS) models were fitted using the GloREDa measured values as an independent variable (Fig. 4). An optimal model should have a zero valued intercept and a regression coefficient as close as possible to one, thus matching the grey line of Fig. 4. Given these criteria our study outperforms the other models³⁰⁻³² with a regression line very close to a coefficient of one. All the other models clearly show a significant bias (over-estimating the predictions) in either the slope or the intercept or both. Moreover, all the other models suffer from high variability as evidenced by the dispersion of the points cloud (Fig. 4).

Table 1 shows the comparison of the fitted models when the full measured range is considered. The intercept of the four models is shown in the first line as “B” column while the second line show the regression coefficients for each of the global models. The past studies³⁰⁻³² have a high intercept bias with values deviating from 0 by more than 1,000, whereas the model proposed in this study has a much smaller value for the intercept (−204.2). The regression coefficient of the current study has a relatively small deviation from 1 (0.13), while the other models have a deviation between −0.27 and −0.85. Also in terms of R² our model clearly outperforms all other models as it was fitted directly to the measured data (Table 1). Solely, the model from Nachtergaele et al.³⁰, performs remarkably well compared to the models of Yang et al.³² and Naipal et al.³¹.

**Sources of uncertainty.** Our inclusive approach to compile the maximum possible number of erosive events assumes that the data collected in one period is comparable with data from other periods. The inter-comparison of different time periods and the non-existence of other alternatives has been extensively discussed²⁸, and has also been followed in similar data collection exercises such as the updated world map of the Koppen-Geiger climate classifications²⁹. Obviously, this simplifying assumption is not valid at local scale, where observed trends in erosivity have been measured over a considerable period (e.g. 30 years). However, at the global scale, the inclusion of 3,625 stations smooths the potential bias due to the presence of trends and long-term temporal variability. Further, we preferred this inclusive approach over the alternative of selecting a common measurement
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of predicted vs. measured R-factor values (values below 10,000 MJ mm⁻¹ ha⁻¹ yr⁻¹) for the three previous and the presented global models. Grey line is the result of an optimal model (Intercept = 0 and regression coefficient = 1); Blue line is the regression result of each model; Grey zone is the 99% confidence interval for the coefficient.</figcaption></figure>

period, since the latter would significantly reduce the station density entailing a reduction of the spatial prediction accuracy.

Past studies recommended a minimum period of 22-years for calculating long term R-factor, while 10-years measurements may lead to under or over-estimation<sup>25, 34, 35</sup>. As 12% of GloREDa includes stations with less than 10-years measurements, we recognize that this may cause a bias due to the temporal variation of R-factor. This bias, even if we are convinced that it is limited, should be considered.

Other uncertainties are related to the methodological approach such as the a) discrepancies in case of using a different algorithm for rainfall energy estimation (\>97.7% of the calculated R-factor stations are based on the original RUSLE equation); b) the application of temporal resolution calibration factors found in European Union to the global dataset; c) the under-representativeness of measured R-factor data in highlands (11.6% of the total dataset is located in areas > 1000 m a.s.l. while the respective area amounts up to 20%), tropical and arid areas and d) the Gaussian Process Regression interpolation model. With regard to the first point the algorithm for unit rainfall energy (Equation 2) may underestimate the rainfall erosivity for instance in high erosive regions (e.g. Ethiopian highlands) due to large drop sizes<sup>36</sup>. The effect of the above listed uncertainties is considered not significant for the objectives of this product which is the model application at global scale. However, we propose that further improvements can be implemented in regional studies taking into account the best suitable algorithm for rainfall unit energy estimation, using more regional-based calibration factors to account for time resolution discrepancies and having a more representative (in terms of climate and topography) pool of stations with measured R-factor data.

**Implications of global erosivity map and data availability.** The new global erosivity map is proposed for global and continental assessments of soil erosion by water, flood risk and natural hazard prevention. Therefore, the aim of the global erosivity dataset is not to challenge other local (or national) erosivity maps, developed from local data with better quality that may not have been available for the present study. Nevertheless, our global R-factor map can potentially cover gaps, where erosivity has not been estimated (due to lack of data), or where it has been calculated solely from rainfall amounts.

Current global estimates on soil erosion by water are very uncertain ranging over one order of magnitude from ca. 20 to over 200 Pg yr⁻¹. More accurate global predictions of the rill and interrill soil erosion rates can only be achieved when the RUSLE rainfall erosivity factor is thoroughly computed. In this study, we present a robust global rainfall-runoff erosivity factor to measure the erosive potential of rainfall, which is a key input in soil erosion prediction models. The global erosivity map represents the first assessment that is solely based on the original
<table><caption>**Table 1.** Ordinary least squares coefficients (B) of global models for the assessment of rainfall erosivity (R-factor).</caption><thead><tr><th></th><th colspan="3">Naipal et al.<sup>31</sup></th><th colspan="3">Nachtergaele et al. (2010)</th><th colspan="3">Yang et al.<sup>32</sup></th><th colspan="3">Panagos et al.</th></tr><tr><th></th><th>B</th><th>std. Error</th><th>p-value</th><th>B</th><th>std. Error</th><th>p-value</th><th>B</th><th>std. Error</th><th>p-value</th><th>B</th><th>std. Error</th><th>p-value</th></tr></thead><tbody><tr><td>Intercept</td><td><strong>1458.2</strong></td><td>279.3</td><td>&lt;0.001</td><td><strong>-1149.0</strong></td><td>141.0</td><td>&lt;0.001</td><td><strong>1223.8</strong></td><td>225.7</td><td>&lt;0.001</td><td><strong>-204.2</strong></td><td>81.2</td><td>0.012</td></tr><tr><td>Regression coefficient</td><td><strong>0.15</strong></td><td>0.12</td><td>0.187</td><td><strong>0.73</strong></td><td>0.04</td><td>&lt;0.001</td><td><strong>0.20</strong></td><td>0.08</td><td>.009</td><td><strong>1.13</strong></td><td>0.05</td><td>&lt;0.001</td></tr><tr><td>Observations</td><td colspan="3">3530</td><td colspan="3">3530</td><td colspan="3">3530</td><td colspan="3">3530</td></tr><tr><td>R<sup>2</sup></td><td colspan="3"><strong>0.155</strong></td><td colspan="3"><strong>0.385</strong></td><td colspan="3"><strong>0.207</strong></td><td colspan="3"><strong>0.811</strong></td></tr></tbody></table>

RUSLE approach using sub-hourly measured rainfall data for 3,625 stations, providing a solid and harmonised basis for a robust spatial interpolation results. Our results show new insights into the global geography of rainfall erosivity and the global erosivity map which will be publicly available can be employed by other research groups to perform national, continental and global soil erosion modelling.

The GloREDa can be considered as an important step to bring together a large group of scientists to advance your understanding of large scale patterns related to land degradation, facing the current policy challenges and demands by the Food and Agriculture Organization (FAO), the Intergovernmental Platform on Biodiversity and Ecosystem Services (IPBES), and the United Nations Convention to Combat Desertification (UNCCD).

The Global erosivity map (GeoTIFF format) at 30 arc-seconds (~1 km) resolution is available for free download in the European Soil Data Centre (ESDAC) website at http://esdac.jrc.ec.europa.eu. The calculated erosivity values per station in GloREDa will become available in the future pending on the agreed copyright issues with data providers.

## Methods

The generation of the global erosivity map involved the following steps: a) the collection of high-temporal resolution rainfall data; b) the calculation of the erosivity factor (R-factor) for each rainfall station; c) the normalisation of R-factor values calculated with rainfall data collected at different time steps (1 min to 60 min), and d) the spatial interpolation of the R-factor point values.

**Collection of data.** At global scale, this is the first time a data collection of observed (measured) high temporal resolution rainfall data takes place. The collection of high temporal resolution rainfall data from the maximum possible number of countries was considered necessary in order to have a representative sample across climatic and geographic gradients.

The Global Rainfall Erosivity Database (GloREDa) was compiled based on the following three components:

* Rainfall Erosivity database at European Scale (REDES) which includes 1,675 rainfall stations in the European Union and Switzerland (28 countries in total). REDES has been the data source for the erosivity map of Europe²⁸ and was updated in 2015 to estimate monthly erosivity in Europe¹³, ³³.

* A global data collection of both high resolution rainfall data and calculated erosivity values based on high temporal resolution rainfall data (Table 2). The data collection yielded 1,865 additional rainfall stations (around 52% of the total stations in GloREDa) from 23 countries outside of Europe. For the majority of the stations (located in China, Japan, India, Kuwait, Israel, Turkey, the United States of America, Mexico, Costa Rica, Jamaica, South Africa and Suriname) the calculation of the rainfall erosivity, based on high temporal resolution data, was performed for the first time (the calculation of rainfall erosivity, named as R-factor, is presented in the homonymous sub-section). For other regions, we used published erosivity data that were based on high temporal resolution rainfall datasets. This included countries such as Australia³⁷, New Zealand²⁷, South Korea³⁸, Iran³⁹, Malaysia⁴⁰, Colombia²¹, Brazil²², Chile²³, Mauritius⁴¹ and Algeria⁴².

* A literature review was used to fill some important data gaps, mainly in Africa, where high temporal resolution rainfall data are scarce. As a result of this exercise, rainfall erosivity values for 85 stations (2.4% of the whole database) from 13 countries were inserted in GloREDa. These countries included Canada⁴³, Argentina²⁴, Cuba⁴⁴, Cape Verde⁴⁵, Cameroon⁴⁶, Eritrea⁴⁷, Ethiopia⁴⁸, Kenya⁴⁹, ⁵⁰, Niger⁵¹, Nigeria⁵², ⁵³, Rwanda⁵⁴, Tenerife⁵⁵ and Zambia⁵⁶.

Summarizing, we collected high temporal resolution rainfall data for 3,540 stations (97.7% of the total GloREDa). The total number of observational years equals 59,380; this results in an average of 16.8 years of high temporal resolution rainfall data per station. As such, GloREDa is the most comprehensive global database including the largest possible number of stations with high temporal resolution rainfall data.

**Calculation of rainfall erosivity (R-factor).** Rainfall erosivity accounts for the combined effect of rainfall duration, magnitude and intensity. In addition, it is also necessary to take into account the frequency of erosive events over a longer time period. In this study, the original R-factor from the Revised Universal Soil Loss Equation (RUSLE)²⁵ was used for the vast majority (>97.7%) of the rainfall stations included in GloREDa. Accordingly, the calculation of rainfall erosivity (EI₃₀) of a single event was based on the following equation:
<table><thead><tr><th>Country</th><th></th><th>No. of Stations</th><th>(Main) Period Covered</th><th>Years per station (average)</th><th>(Main) Temporal resolution of rainfall data</th><th>Source of high temporal resolution rainfall data</th></tr></thead><tbody><tr><td>AT</td><td>Austria</td><td>84</td><td>1995–2010</td><td>21</td><td>27 stations: 10Min; 57 stations: 15Min</td><td>Hydrographic offices of Upper Austria, Lower Austria, Burgenland, Styria, Salzburg, Carinthia, Vorarlberg and Tyrol.</td></tr><tr><td>AU</td><td>Australia</td><td>167</td><td>1961–2000</td><td>40</td><td>6 Min</td><td>Bureau of Meteorology in Australia; Yu et al.<sup>37</sup></td></tr><tr><td>BE</td><td>Belgium - Flanders</td><td>20</td><td>2004–2013</td><td>10</td><td>30 Min</td><td>Flemish Environmental Agency (VMM),</td></tr><tr><td>BE</td><td>Belgium - Wallonia</td><td>29</td><td>2004–2013</td><td>10</td><td>60 Min</td><td>Service Public de Wallonie</td></tr><tr><td>BG</td><td>Bulgaria</td><td>84</td><td>1951–1976</td><td>26</td><td>30 Min</td><td>Rousseva et al.<sup>37</sup></td></tr><tr><td>BR</td><td>Brazil</td><td>87</td><td>1986–2008</td><td>19</td><td>49 stations: 5 Min; 38 stations: 10 Min</td><td>Oliveira et al.<sup>22</sup> and updated in 05/2016</td></tr><tr><td>CH</td><td>Switzerland</td><td>71</td><td>1988–2010</td><td>22</td><td>10 Min</td><td>MeteoSwiss</td></tr><tr><td>CL</td><td>Chile</td><td>18</td><td>1976–1995</td><td>17</td><td>15 Min</td><td>General Directorate of Water Resources (DGA), Government of Chile</td></tr><tr><td>CN</td><td>China</td><td>387</td><td>1998–2012</td><td>14</td><td>60 Min</td><td>National Meteorological Information Center of the China Meteorological Administration</td></tr><tr><td>CO</td><td>Colombia</td><td>6</td><td>1987–1996</td><td>10</td><td>5 Min</td><td>Centro Nacional de Investigaciones de Café - Cenicafé</td></tr><tr><td>CR</td><td>Costa Rica</td><td>5</td><td>2011–2015</td><td>6</td><td>30 Min</td><td>University of Costa Rica, Costa Rica</td></tr><tr><td>CY</td><td>Cyprus</td><td>35</td><td>1974–2013</td><td>39</td><td>30 Min</td><td>Cyprus Department of Meteorology</td></tr><tr><td>CZ</td><td>Czech Republic</td><td>32</td><td>1961–1999</td><td>35</td><td>30 Min</td><td>Research Institute for Soil and Water Conservation (Czech Republic)</td></tr><tr><td>DE</td><td>Germany</td><td>148</td><td>1996–2013</td><td>18</td><td>60 Min</td><td>Deutscher Wetterdienst (DWD)</td></tr><tr><td>DK</td><td>Denmark</td><td>30</td><td>1988–2012</td><td>15</td><td>60 Min</td><td>Danish Meteorological Institute (DMI), Aarhus University</td></tr><tr><td>DZ</td><td>Algeria</td><td>120</td><td>1977–2004</td><td>24</td><td>15 Min</td><td>National Agency of Hydraulic Resources, Algeria</td></tr><tr><td>EE</td><td>Estonia</td><td>21</td><td>2007–2013</td><td>7</td><td>60 Min</td><td>Estonian Environment Agency</td></tr><tr><td>ES</td><td>Spain</td><td>146</td><td>2002–2013</td><td>12</td><td>24 stations: 10 Min; 104 stations: 15 Min; 18 stations 30 MIN</td><td>Regional water agencies</td></tr><tr><td>FI</td><td>Finland</td><td>64</td><td>2007–2013</td><td>7</td><td>60 Min</td><td>Finnish Climate Service Centre (FMI)</td></tr><tr><td>FR</td><td>France</td><td>81</td><td>2004–2013</td><td>10</td><td>60 Min</td><td>Météo-France DP/SERV/FDP</td></tr><tr><td>GR</td><td>Greece</td><td>80</td><td>1974–1997</td><td>30</td><td>30 Min</td><td>Hydroskopio</td></tr><tr><td>HR</td><td>Croatia</td><td>42</td><td>1961–2012</td><td>40</td><td>10 Min</td><td>Croatian Meteo & Hydrological Service</td></tr><tr><td>HU</td><td>Hungary</td><td>30</td><td>1998–2013</td><td>16</td><td>10 Min</td><td>Hungarian Meteorological Service</td></tr><tr><td>IE</td><td>Ireland</td><td>13</td><td>1950–2010</td><td>56</td><td>60 Min</td><td>Met Éireann - The Irish National Meteorological Service</td></tr><tr><td>IL</td><td>Israel</td><td>61</td><td>1998–2015</td><td>17</td><td>30 Min</td><td>Israel Meteorological Service</td></tr><tr><td>IN</td><td>India</td><td>247</td><td>2007–2015</td><td>7</td><td>60 Min</td><td>India Meteorological Department, Ministry of Earth Sciences</td></tr><tr><td>IR</td><td>Iran</td><td>70</td><td>1984–2004</td><td>21</td><td>10 Min</td><td>Iranian Meteorological Organization</td></tr><tr><td>IT</td><td>Italy</td><td>251</td><td>2002–2011</td><td>10</td><td>30 Min</td><td>Regional meteorological services, Regional agencies for environmental protection (ARPA)</td></tr><tr><td>JM</td><td>Jamaica</td><td>9</td><td>2003–2014</td><td>12</td><td>2 Min</td><td>Meteorological service Jamaica</td></tr><tr><td>JP</td><td>Japan</td><td>55</td><td>2006–2015</td><td>10</td><td>60 Min</td><td>Japan Meteorological Agency (JMA)</td></tr><tr><td>KR</td><td>South Korea</td><td>75</td><td>1998–2015</td><td>18</td><td>10 Min</td><td>Korea Meteorological Administration (KMA)</td></tr><tr><td>KW</td><td>Kuwait</td><td>15</td><td>2007–2015</td><td>9</td><td>60 Min</td><td>Department of Meteorology, Directorate General of Civil Aviation, State of Kuwait</td></tr><tr><td>LT</td><td>Lithuania</td><td>3</td><td>1992–2007</td><td>16</td><td>30 Min</td><td>Lithuanian Agriculture and Forestry Research Centre</td></tr><tr><td>LU</td><td>Luxembourg</td><td>16</td><td>2000–2013</td><td>11</td><td>60 Min</td><td>Agrarmeteorologisches Messnetz</td></tr><tr><td>LV</td><td>Latvia</td><td>4</td><td>2007–2013</td><td>7</td><td>60 Min</td><td>Latvian Environment, Geology and Meteorology Centre</td></tr><tr><td>MU</td><td>Mauritius</td><td>5</td><td>2005–2008</td><td>5</td><td>6 Min</td><td>Mauritius Meteorological Services (MMS)</td></tr><tr><td>MX</td><td>Mexico</td><td>15</td><td>2006–2014</td><td>9</td><td>60 Min</td><td>CONAGUA, Comisión Nacional Del Agua, Servicio Meteorológico Nacional, Mexico.</td></tr><tr><td>MY</td><td>Malaysia</td><td>2</td><td>1981–1998</td><td>18</td><td>10 Min</td><td>Yu et al.<sup>40</sup></td></tr><tr><td>NL</td><td>Netherlands</td><td>32</td><td>1981–2010</td><td>24</td><td>60 Min</td><td>Royal Netherlands Meteorological Institute</td></tr></tbody><tfoot><tr><td colspan="7">Continued</td></tr></tfoot></table>
<table><caption>**Table 2.** Overview of the high resolution rainfall data used to estimate global rainfall erosivity. In addition, erosivity information of 85 stations from 13 countries found in the literature<sup>24, 43–56</sup> was included in the global map (not shown in the table).</caption><thead><tr><th colspan="2">Country</th><th>No. of Stations</th><th>(Main) Period Covered</th><th>Years per station (average)</th><th>(Main) Temporal resolution of rainfall data</th><th>Source of high temporal resolution rainfall data</th></tr></thead><tbody><tr><td>NZ</td><td>New Zealand</td><td>35</td><td>2000–2012</td><td>12</td><td>10 Min</td><td>New Zealand Institute of Water and Atmospheric Research (NIWA)</td></tr><tr><td>PL</td><td>Poland</td><td>13</td><td>1961–1988</td><td>27</td><td>30 Min</td><td>Warsaw University of Life Sciences</td></tr><tr><td>PT</td><td>Portugal</td><td>41</td><td>2001–2012</td><td>11</td><td>60 Min</td><td>Agência Portuguesa do Ambiente</td></tr><tr><td>RO</td><td>Romania</td><td>60</td><td>2006–2013</td><td>8</td><td>10 Min</td><td>Meteorological Administration</td></tr><tr><td>RU</td><td>Russian Federation</td><td>218</td><td>1961–1983</td><td>23</td><td>30 Min</td><td>Lomonosov Moscow State University</td></tr><tr><td>SE</td><td>Sweden</td><td>73</td><td>1996–2013</td><td>18</td><td>60 Min</td><td>Swedish Meteorological and Hydrological Institute (SMHI)</td></tr><tr><td>SI</td><td>Slovenia</td><td>31</td><td>1999–2008</td><td>10</td><td>5 Min</td><td>Slovenian Environment Agency</td></tr><tr><td>SK</td><td>Slovakia</td><td>103</td><td>1971–1990</td><td>20</td><td>60 Min</td><td>Slovak Hydrometeorological Institute, Climatological service</td></tr><tr><td>SR</td><td>Suriname</td><td>11</td><td>1987–2010</td><td>25</td><td>60 Min</td><td>Meteorological organization of Suriname</td></tr><tr><td>TR</td><td>Turkey</td><td>160</td><td>2005–2014</td><td>10</td><td>1 Min</td><td>Ministry of Forest and Water Affairs General Directorate of Combating Desertification and Erosion</td></tr><tr><td>UK</td><td>United Kingdom</td><td>38</td><td>1993–2012</td><td>20</td><td>60 Min</td><td>NERC &amp; UK Environ. Change Network(ECN), British Atmospheric Data Centre (BADC)</td></tr><tr><td>US</td><td>United States of America</td><td>92</td><td>2006–2016</td><td>11</td><td>5 Min</td><td>U.S. Climate Reference Network (USCRN), NOAA; Diamond et al.<sup>58</sup></td></tr><tr><td>ZA</td><td>South Africa</td><td>5</td><td>2001–2005</td><td>5</td><td>5 Min</td><td>Nel and Summer<sup>41</sup></td></tr><tr><td>Total</td><td></td><td>3,540</td><td></td><td></td><td></td><td></td></tr></tbody></table>

$$EI_{30} = \left( \sum_{r=1}^{k} e_r v_r \right) I_{30} \qquad (1)$$

where $e_r$ is the unit rainfall energy (MJ ha<sup>-1</sup> mm<sup>-1</sup>) and $v_r$ the rainfall volume (mm) during the r<sup>th</sup> time period of a rainfall event divided in k-parts. $I_{30}$ is the maximum 30-minutes rainfall intensity (mm h<sup>-1</sup>). The unit rainfall energy ($e_r$) is calculated for each time interval as follows<sup>59</sup>:

$$e_r = 0.29[1 - 0.72e^{(-0.05i_r)}] \qquad (2)$$

where $i_r$ is the rainfall intensity during the time interval (mm h<sup>-1</sup>).

R is the average annual rainfall erosivity (MJ mm ha<sup>-1</sup> h<sup>-1</sup> yr<sup>-1</sup>):

$$R = \frac{\sum_{j=1}^{n} \sum_{k=1}^{m_j} (EI_{30})_k}{n} \qquad (3)$$

where n is the number of years recorded, $m_j$ is the number of erosive events during a given year j and k is the index of a single event with its corresponding erosivity $EI_{30}$.

Equation (2) was developed by Brown and Foster<sup>59</sup> as part of the Revised Universal Soil Loss Equation (RUSLE, Renard et al.<sup>25</sup>), and replaced the original equation of Wischmeier and Smith<sup>35</sup> used in the Universal Soil Loss Equation (USLE). Equation (2) was further modified as part of RUSLE2<sup>60</sup>, but its use has been mostly limited to the United States. However, outside of the United States the RUSLE equation is still the most commonly used and we have therefore applied it for our erosivity calculations. A comparison of the two equations yields slightly higher erosivity values for the RUSLE2 equation<sup>60</sup>. A recent study showed that empirical rainfall kinetic energy relationships compare well to measurement when complete events were considered (R<sup>2</sup> > 0.90). Nonetheless, future research will explore how our current GloREDa results compare with alternative functions for rainfall energy.

According to the RUSLE handbook<sup>25</sup>, the erosive rainfall events were computed based on the following criteria: (i) the cumulative rainfall of an event is greater than 12.7 mm, (ii) the event has at least one peak that is greater than 6.35 mm during a period of 15 min (or 12.7 mm during a period of 30 min) and, (iii) a rainfall accumulation of less than 1.27 mm during a period of six hours splits a longer storm period into two storms. The erosivity factor equations and the above mentioned criteria have been developed according to more than 10,000 plot-years of experiments. The R-factor for each station in GloREDa was calculated using the Rainfall Intensity Summarisation Tool (RIST) software developed by the United States Department of Agriculture (USDA).

**Calibration of different time resolutions.** For such an extensive data collection, it is expected to have a variety in both i) the range of available data-years and ii) the time resolution of the data. According to the
GloREDa statistics, 35.7% of the stations had rainfall data at very high resolution (<=15 min); ~25.7% had an intermediate resolution (30 minutes), while the remaining 38.6% had a resolution of 60 minutes. Due to this heterogeneity, a time step calibration of erosivity values was considered necessary.

We selected a 30-minute time resolution as an acceptable compromise between the coarse time resolution of 60 minutes and the higher ones (<=15 minutes). For the calibration, we have selected a pool of 86 rainfall stations as their data were available at multiple temporal resolutions and their geographical coverage is representing 14 countries covering a large climate gradient. The calibration process included three steps: a) The R-factor was calculated at the highest possible resolution (5, 10 or 15 minutes); b) The rainfall data were aggregated to coarser resolutions (30 and 60 minutes) and the corresponding R-factor was calculated at the coarser resolutions and c) The calibration factors was derived from the regression analysis of the R-factor calculations at the highest and coarser resolutions. Those calibration factors have been developed in the European study³³ and were in agreement with range values provided in the literature²⁵,⁶⁰,⁶¹ which have been calculated in China, South Italy and U.S.A.

**Interpolation of GloREDa.** Given the correlation between the R-factor and monthly climatic data²⁸, a regression approach was used to infer the spatial global distribution of rainfall erosivity from a series of independent climatic covariates derived from the WorldClim database⁶². The covariates used included average monthly precipitation, average minimum and maximum monthly precipitation, average monthly temperature, precipitation of the wettest month, precipitation of the driest month, and precipitation seasonality, as defined in the WorldClim database (www.worldclim.org). These variables represent long term conditions based on the interpolation of observed data for the period 1960–1990⁶². As a result, our erosivity map represents long-term erosivity patterns. Elevation was not included in the model, as it was already used to estimate some of the WorldClim variables. We subsequently assembled a dataset where each location (i.e. rainfall station) had a mean R-factor (independent variable) as well as values for the climatic variables (independent variables), and used it as input for the interpolation.

We used Gaussian Process Regression (GPR)⁶³,⁶⁴ due to the large number of support covariates, their potential collinearity, and the presence of non-linear relationships between the target variable (R-factor) and the covariates. GPR is a non-linear regression approach that can model non-linear relations by projecting the inputs into a higher dimensional space using basis functions, and by creating a regression model in said space. Considering the regression form $y = f(x) + N(0, \sigma_n^2)$ with $f(x) = x^T w$ (where $y$ are the observed responses, $x$ the covariates values vector, $f$ a set of functions and $w$ a vector of weights), the GPR uses a projection from the original input space into a feature space using kernel expansion so that $f(x)$ can be rewritten as $f(x) = \phi(x)^T w$ where $\phi$ is the kernel function.

In this study the Radial Basis Function (RBF) Gaussian kernel was used, which can be written as
$$K(\mathbf{x}, \mathbf{x}') = \exp\left(-\frac{\|\mathbf{x} - \mathbf{x}'\|^2}{2\sigma^2}\right)$$ where $\sigma$ is free tunable parameter. The Gaussian kernel is highly adaptable and commonly applied in machine learning⁶⁵. The main advantages of GPR are that it can model very complex non-linear relations between covariates and the target variable, and directly model both average and variance estimation, thus providing information about prediction uncertainty. Moreover, GPR is resistant to the issues derived from collinearity among independent variables that can arise in other statistical models. Finally, the GPR has the advantage of adapting in local conditions (climatic) due to its non-parametric nature. In this study the GPR model performance was tested for both, a fitting, and a cross-validation dataset. The cross-validation was carried out by random sampling with 10% replacement of the original dataset used for validation.

**References**

1. Freebairn, J. W. & Zillman, J. W. Economic benefits of meteorological services. *Meteorol. Appl.* **9**, 33–44 (2002).

2. Tapiador, F. J. et al. Global precipitation measurement: Methods, datasets and applications. *Atmospheric Res.* **104**, 70–97 (2012).

3. Alexander, L. V. et al. Global observed changes in daily climate extremes of temperature and precipitation. *J. Geophys. Res. Atmospheres* **111** (2006).

4. Orlowsky, B. & Seneviratne, S. I. Global changes in extreme events: regional and seasonal dimension. *Clim. Change* **110**, 669–696 (2012).

5. Olsson, J., Berggren, K., Olofsson, M. & Viklander, M. Applying climate model precipitation scenarios for urban hydrological assessment: A case study in Kalmar City, Sweden. *7th Int. Workshop Precip. Urban Areas7th Int. Workshop Precip. Urban Areas* **92**, 364–375 (2009).

6. Willems, P., Arnbjerg-Nielsen, K., Olsson, J. & Nguyen, V. T. V. Climate change impact assessment on urban rainfall extremes and urban drainage: methods and shortcomings. *Atmospheric Res.* **103**, 106–118 (2012).

7. Olesen, J. E. et al. Impacts and adaptation of European crop production systems to climate change. *Eur. J. Agron.* **34**, 96–112 (2011).

8. Ramon, R., Minella, J. P. G., Merten, G. H., de Barros, C. A. P. & Canale, T. Kinetic energy estimation by rainfall intensity and its usefulness in predicting hydrosedimentological variables in a small rural catchment in southern Brazil. *Achiev. Sustain. Conserv. Chang. World* **148**, Part 2, 176–184 (2017).

9. Vereecken, H. et al. Modeling Soil Processes: Review, Key Challenges, and New Perspectives. *Vadose Zone J.* **15** (2016).

10. UNCCD. *United Nations Convention to Combat Desertification: Zero Net Land Degradation*. (Bonn, 2012).

11. IPBES. *First Order Drafts of IPBES assessments*. (IPBES, 2016).

12. FAO and ITPS. *Status of the World's Soil Resources (SWSR) – Main Report*. (Food and Agriculture Organization of the United Nations and Intergovernmental Technical Panel on Soils, 2015).

13. Ballabio, C. et al. Mapping monthly rainfall erosivity in Europe. *Sci. Total Environ.* **579**, 1298–1315 (2017).

14. Zipser, E. J., Cecil, D. J., Liu, C., Nesbitt, S. W. & Yorty, D. P. Where are the most intense thunderstorms on Earth? *Bull. Am. Meteorol. Soc.* **87**, 1057 (2006).

15. Diodato, N., Knight, J. & Bellocchi, G. Reduced complexity model for assessing patterns of rainfall erosivity in Africa. *Glob. Planet. Change* **100**, 183–193 (2013).

16. Vrieling, A., Hoedjes, J. C. & van der Velde, M. Towards large-scale monitoring of soil erosion in Africa: Accounting for the dynamics of rainfall erosivity. *Glob. Planet. Change* **115**, 33–43 (2014).

17. Zhu, Z. & Yu, B. Validation of rainfall erosivity estimators for mainland China. *Trans. ASABE* **58**, 61–71 (2015).
18. Qin, W. et al. Spatial distribution and temporal trends of rainfall erosivity in mainland China for 1951–2010. *Catena* **147**, 177–186 (2016).

19. Shiono, T., Ogawa, S., Miyamoto, T. & Kameyama, K. Expected impacts of climate change on rainfall erosivity of farmlands in Japan. *Ecol. Eng.* **61**, 678–689 (2013).

20. Leow, C. S., Ghani, A. B., Zakaria, N. A. & Abidin, R. Z. Development of rainfall erosivity isohyet map for Peninsular Malaysia. in (2011).

21. Hoyos, N., Waylen, P. R. & Jaramillo, Á. Seasonal and spatial patterns of erosivity in a tropical watershed of the Colombian Andes. *J. Hydrol.* **314**, 177–191 (2005).

22. Oliveira, P. T. S., Wendland, E. & Nearing, M. A. Rainfall erosivity in Brazil: A review. *CATENA* **100**, 139–147 (2013).

23. Bonilla, C. A. & Vidal, K. L. Rainfall erosivity in Central Chile. *J. Hydrol.* **410**, 126–133 (2011).

24. Rojas, A. & Conde, A. Estimación del factor R de la ecuación universal de pérdidas de suelo para el centro-este de la República Argentina. *Cienc. Suelo* **3**, 85–94 (1985).

25. Renard, K. G., Foster, G. R., Weesies, G. A., McCool, D. K. & Yoder, D. C. Predicting soil erosion by water: a guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). **703**, (US Government Printing Office Washington, DC, 1997).

26. Teng, H. et al. Assimilating satellite imagery and visible-near infrared spectroscopy to model and map soil loss by water erosion in Australia. *Environ. Model. Softw.* **77**, 156–167 (2016).

27. Klik, A., Haas, K., Dvorackova, A. & Fuller, I. C. Spatial and temporal distribution of rainfall erosivity in New Zealand. *Soil Res.* **53**, 815–825 (2015).

28. Panagos, P. et al. Rainfall erosivity in Europe. *Sci. Total Environ.* **511**, 801–814 (2015).

29. Peel, M. C., Finlayson, B. L. & McMahon, T. A. Updated world map of the Köppen-Geiger climate classification. *Hydrol. Earth Syst. Sci. Discuss.* **4**, 439–473 (2007).

30. Nachtergaele, F. O. et al. Global Land Degradation Information System (GLADIS), an information database for land degradation assessment at global level. (FAO, 2011).

31. Naipal, V., Reick, C., Pongratz, J. & Van Oost, K. Improving the global applicability of the RUSLE model – adjustment of the topographical and rainfall erosivity factors. *Geosci Model Dev* **8**, 2893–2913 (2015).

32. Yang, D., Kanae, S., Oki, T., Koike, T. & Musiak, K. Global potential soil erosion with reference to land use and climate changes. *Hydrol. Process.* **17**, 2913–2928 (2003).

33. Panagos, P. et al. Monthly rainfall erosivity: conversion factors for different time resolutions and regional assessments. *Water* **8**, 119 (2016).

34. Verstraeten, G., Poesen, J., Demarée, G. & Salles, C. Long-term (105 years) variability in rain erosivity as derived from 10-min rainfall depth data for Ukkel (Brussels, Belgium): Implications for assessing soil erosion rates. *J. Geophys. Res. Atmospheres* **111**, n/a-n/a (2006).

35. Wischmeier, W. H. & Smith, D. D. Predicting rainfall erosion losses: A guide to conservation planning. *Predict. Rainfall Eros. Losses-Guide Conserv. Plan.* (1978).

36. Nyssen, J. et al. Rainfall erosivity and variability in the Northern Ethiopian Highlands. *J. Hydrol.* **311**, 172–187 (2005).

37. Yu, B. An assessment of uncalibrated CLIGEN in Australia. *Agric. For. Meteorol.* **119**, 131–148 (2003).

38. Risal, A. et al. Application of Web ERosivity Module (WERM) for estimation of annual and monthly R factor in Korea. *CATENA* **147**, 225–237 (2016).

39. Sadeghi, S. H. R. & Hazbavi, Z. Trend analysis of the rainfall erosivity index at different time scales in Iran. *Nat. Hazards* **77**, 383–404 (2015).

40. Yu, B., Hashim, G. M. & Eusof, Z. Estimating the R-factor with limited rainfall data: a case study from peninsular Malaysia. *J. Soil Water Conserv.* **56**, 101–105 (2001).

41. Nel, W. et al. Intra-event characteristics of extreme erosive rainfall on Mauritius. *Phys. Geogr.* **37**, 264–275 (2016).

42. Meddi, M., Toumi, S. & Assani, A. A. Spatial and temporal variability of the rainfall erosivity factor in Northern Algeria. *Arab. J. Geosci.* **9**, 282 (2016).

43. Wall, G. J., Coote, D. R., Pringle, E. A. & Shelton, I. J. Revised Universal Soil Loss Equation for Application in Canada: A Handbook for Estimating Soil Loss from Water Erosion in Canada. *Agric. Agri-Food Can. Res. Branch Ott. Contrib. No AAFAAC2244E* (2002).

44. Schiettecatte, W. et al. Influence of landuse on soil erosion risk in the Cuyaguateje watershed (Cuba). *Catena* **74**, 1–12 (2008).

45. Mannnaerts, C. & Gabriels, D. Rainfall erosivity in Cape Verde. *Soil Tillage Res.* **55**, 207–212 (2000).

46. Fournier, J. Agressivité climatique et risques érosifs dans la région de Dschang, ouest Cameroun. *Bull. Réseau Eros.* **14**, 145–156 (1994).

47. Stillhardt, B., Herweg, K. & Hurni, H. Long-term Monitoring of Soil Erosion and Soil and Water Conservation in Afdeyu, Eritrea (1984–1998). *CDE Univ. Berne Berne Switz.* (2002).

48. Haile, M., Herweg, K. & Stillhardt, B. *Sustainable Land Management – A New Approach to Soil and Water Conservation in Ethiopia*. (Mekelle University (with CDE-University of Berne and NCCR), 2006).

49. Angima, S., Stott, D., O'Neill, M., Ong, C. & Weesies, G. Soil erosion prediction using RUSLE for central Kenyan highland conditions. *Agric. Ecosyst. Environ.* **97**, 295–308 (2003).

50. Ulsaker, L. G. & Onstad, C. A. Relating rainfall erosivity factors to soil loss in Kenya. *Soil Sci. Soc. Am. J.* **48**, 891–896 (1984).

51. Delwaulle, J. C. Résultats de six ans d'observations sur l'érosion au Niger. *Bois For. Trop.* **150**, 15–36 (1973).

52. Gobin, A. M., Campling, P., Deckers, J. A., Poesen, J. & Feyen, J. Soil erosion assessment at the Udi-Nsukka Cuesta (southeastern Nigeria). *Land Degrad. Dev.* **10**, 141–160 (1999).

53. Salako, F. K., Ghuman, B. S. & Lal, R. Rainfall erosivity in south-central Nigeria. *Soil Technol.* **7**, 279–290 (1995).

54. Ryumugabe, J. B. & Berding, F. R. Variabilité de l'indice d'agressivité des pluies au Rwanda. *Bull. Réseau Eros.* **12**, 113–119 (1992).

55. Rodriguez, A. et al. Forms of eroded soil organic carbon in andosols of the Canary Islands (Spain). *Geoderma* **121**, 205–219 (2004).

56. Pauwelyn, P. L. L., Lenvain, J. S. & Sakala, W. K. Iso-erodent map of Zambia. *Soil Technol.* **1**, 235–250 (1988).

57. Rousseva, S. et al. Soil erosion risk in Bulgaria and recommendations for soil protective use of agricultural lands. *Part II South. Bulg. Agric. Acad. Inst. Soil Sci. Nikola Poushkarov Publ.-Eco Sofia Bulg.* (2010).

58. Diamond, H. J. et al. U.S. Climate Reference Network after One Decade of Operations: Status and Assessment. *Bull. Am. Meteorol. Soc.* **94**, 485–498 (2013).

59. Brown, L. C. & Foster, G. R. Storm Erosivity Using Idealized Intensity Distributions. *Trans. Am. Soc. Agric. Eng.* **30**, 379–386 (1987).

60. Yin, S., Xie, Y., Liu, B. & Nearing, M. A. Rainfall erosivity estimation based on rainfall data collected over a range of temporal resolutions. *Hydrol. Earth Syst. Sci.* **19**, 4113–4126 (2015).

61. Porto, P. Exploring the effect of different time resolutions to calculate the rainfall erosivity factor R in Calabria, southern Italy. *Hydrol. Process.* **30**, 1551–1562 (2016).

62. Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. *Int. J. Climatol.* **25**, 1965–1978 (2005).

63. Rasmussen, C. E. & Williams, C. K. I. Gaussian processes for machine learning. (MIT Press, 2006).

64. Stein, M. L. Interpolation of spatial data: some theory for kriging. (Springer Science & Business Media, 2012).

65. Hofmann, T., Schölkopf, B. & Smola, A. J. Kernel methods in machine learning. *Ann. Stat.* **117**–1220 (2008).
## Acknowledgements

The authors would also like to acknowledge the following services for proving their data: Bureau of Meteorology in Australia, New Zealand Institute of Water and Atmospheric Research (NIWA), Japan Meteorological Agency (JMA), Korea Meteorological Administration (KMA), National Meteorological Information Center of China, India Meteorological Department, Ministry of Earth Sciences, Iranian Meteorological Organization, Department of Meteorology, Directorate General of Civil Aviation - State of Kuwait, Lomonosov Moscow State University, Israel Meteorological Service, Turkish Ministry of Forestry and Water Affairs, U.S. Climate Reference Network (USCRN, NOAA), Comision Nacional Del Agua, Servicio Meteorologico Nacional, Mexico, Meteorological service Jamaica, University of Costa Rica (UCR), Centro Nacional de Investigaciones de Café - Cenicafé (Colombia), General Directorate of Water Resources (Chile), Meteorological organization of Suriname, Mauritius Meteorological Services (MMS), Algerian National Agency of Hydraulic Resources. Plus all the European Organizations listed in Panagos et al²⁸.

## Author Contributions

Panos Panagos, Cristiano Ballabio, Pasquale Borrelli, Katrin Meusburger elaborated the results, developed the erosivity model and wrote major part of the manuscript. Prof. Bofu Yu, Dr. Andreas Klik, Dr. Kyoung Jae Lim, Dr. Jae E. Yang, Prof. Jinren Ni, Prof. Chiyuan Miao, Dr. Nabansu Chattopadhyay, Prof. Seyed Hamidreza Sadeghi, Dr. Zeinab Hazbavi, Dr. Mohsen Zabihi, Prof. Genady Larionov, Dr. Andrey Gorobets, Dr. Sergey Krasnov, Mr. Yoav Levi, Dr. Gunay Erpul, Prof. Christian Birkel, Dr. Natalia Hoyos, Dr. Victoria Naipal, Dr. Paulo Tarso Oliveira, Prof. Carlos Bonilla, Dr. Mohamed Meddi, Dr. Werner Nel, Mr. Hassan Al Dashti, Mr. Martino Boni, Dr. Nazzareno Diodato, Prof. Kristof Van Oost and Dr. Mark A. Nearing provided the data input, contributed with modelling and analysis and reviewed the manuscript.

## Additional Information

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_6.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017