Rainfall erosivity in Italy: a national scale spatio-temporal assessment

ABSTRACT Soil erosion by water is a serious threat for the Mediterranean region. Raindrop impacts and consequent runoff generation are the main driving forces of this geomorphic process of soil degradation. The potential ability for rainfall to cause soil loss is expressed as rainfall erosivity, a key parameter required by most soil loss prediction models. In Italy, rainfall erosivity measurements are limited to few locations, preventing researchers from effectively assessing the geography and magnitude of soil loss across the country. The objectives of this study were to investigate the spatio-temporal distribution of rainfall erosivity in Italy and to develop a national-scale grid-based map of rainfall erosivity. Thus, annual rainfall erosivity values were measured and subsequently interpolated using a geostatistical approach. Time series of pluviographic records (10-years) with high temporal resolution (mostly 30-min) for 386 meteorological stations were analysed. Regression-kriging was used to interpolate rainfall erosivity values of the meteorological stations to an Italian rainfall erosivity map (500-m). A set of 23 environmental covariates was tested, of which seven covariates were selected based on a stepwise approach (mostly significant at the 0.01 level). The interpolation method showed a good performance for both the cross-validation data set ( = 0.777) and the fitting data set (R2 = 0.779).


Introduction
Soil erosion by water is the ongoing geomorphic process of soil degradation studied most in the Mediterranean region (Kosmas et al. 1997;Boardman and Poesen 2006). It is a widespread phenomenon (Cerdan et al. 2010) that affects several Mediterranean landscapes ) causing on-site and off-site environmental impacts (Rodolfi 2006;Morgan 2009). Recognized in the EU Thematic Strategy for Soil Protection as one of the major threats to European soils (European Commission 2006), soil erosion is a major concern in landscape management and conservation planning (CAP 2003;GAEC 2013).
Raindrop impacts and the consequent runoff generation represent the driving forces of the process (Moss and Green 1983). Soil particle movement occurs when these forces exceed the cohesive forces between the soil particles (Bryan 2000). The potential ability for rainfall to cause soil loss is expressed as rainfall erosivity (Ekern 1954;Wischmeier 1959). The comprehensive rainfall-runoff erosivity factor R, proposed in the USLE and RUSLE soil erosion prediction models (Wischmeier and Smith 1978;Renard et al. 1997), is recognized as one of the best indicators to measure the erosive potential of a rainstorm (Renard et al. 1991). Additionally, the R-factor is also found to be a suitable indicator for hazardous environmental and hydrological processes, for example, landslides and floods (Diodato 2004;Capolongo et al. 2008;Diodato et al. 2015). The R-factor is expressed as the product of the total kinetic energy of a storm times its 30-min maximum rainfall intensity (KE × I 30 ) (Wischmeier and Smith 1978). The best method to calculate the R-factor requires continuous record of subhourly rainfall intensity data (1-5-min intervals) (Williams and Sheridan 1991). However, due to instrumental limitations the R-factor is communally computed at 15-to 30-min intervals. Afterwards conversion factors are generally used to avoid possible underestimations of the R-factor (Williams and Sheridan 1991;Renard et al. 1997;Agnese et al. 2006;Yin et al. 2007). During the last two decades, USLE-based models have become the most applied tools to quantify soil loss rates at multiple scales over Europe. The lack of pluviographic records and the time-consuming procedures required to quantify the R-factor (Diodato 2004) have limited the computation of the USLE rainfall erosivity factor to time-series analyses of single stations in several studies (Zanchi 1988;Verstraeten et al. 2006;Borrelli et al. 2013). Alternatively, for regional-scale studies, approximation equations based on the rainfall volume, instead of the intensity, were adopted for R-factor estimation (Renard and Freimund 1994;De Santos Loureiro and De Azevedo Coutinho 2001;Diodato and Bellocchi 2007;Capolongo et al. 2008;Bosco et al. 2015, among others). In the Algarve region of Portugal, Goovaerts (1999) tested a geostatistical approach based on the incorporation of exhaustive secondary information (elevation) into the spatial prediction of rainfall erosivity. An approach that was further studied and developed with good results by combining rainfall erosivity values measured from high-temporal-resolution records (15-and 30-min) with a set of exhaustive secondary environment variables, for example, the Ebro catchment in Spain (Angulo-Martínez et al. 2009), Switzerland (Meusburger et al. 2012) and Greece (Panagos et al. 2016). However, for large-scale studies (Van der Knijff, Jones, and Montanarella 1999Montanarella , 2000Grimm et al. 2003;Bosco et al. 2015, among others), a simple algorithm assuming that the R-factor is proportional to the cumulative yearly rainfall amount (Diodato and Bellocchi 2007) is generally preferred. This is a simplistic solution based on the rainfall volume, reflecting specific local-type conditions (Goldman, Jackson, and Bursztynsky 1986;Diodato andBellocchi 2007, 2010), which may return only mere approximations of rainfall erosivity.
Studies conducted along the Italian Peninsula indicate that most, if not all, agricultural cropland areas experience some degree of soil erosion by water (APAT 2009). However, measurements of the rainfall erosivity are spatially limited to a few locations, preventing researchers from effectively assessing the geography and magnitude of soil loss across the country. The objectives of this study were to investigate the spatio-temporal distribution of rainfall erosivity in Italy and to develop a nationalscale grid-based map of rainfall erosivity. The procedure consisted of two steps: first, rainfall erosivity was measured following the methods proposed for the (Revised) Universal Soil Loss Equation handbook (Renard et al. 1997). A time series of pluviographic records (10 years) with high temporal resolution (mostly 30-min) for 386 meteorological stations was analysed. Second, a regression-kriging approach was used to interpolate the rainfall erosivity values of the meteorological stations to an Italian rainfall erosivity map at 500-m resolution.
This in-depth study outlines the spatio-temporal patterns of rainfall erosivity across Italy using exclusively measured data, advanced spatial-interpolation techniques and GIS operations. It constitutes a meaningful contribution to Digital Earth, as it explores the use of geo-referenced data to represent the geography of a complex natural phenomenon such as rainfall erosivity at national level. Moreover, the outcomes of this study provide innovative knowledge that can enable more effective assessments of soil loss (multi-scale and multi-temporal).

Data collection
High-temporal-resolution pluviographic records were collected from 386 electronic rain-gauges of the Regional Hydrographic Services and Regional Agencies for Environmental Protection (ARPA). The selected stations are well distributed throughout Italy, covering the different biogeographical zones (EEA 2014) and climate regions (Brunetti et al. 2004). The time-series were composed of continuous 10-year records for most of the stations (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011), that is, the minimum time period to sufficiently represent the inter-annual variability of rainfall erosivity (Renard and Freimund 1994;da Silva 2004).
Sixteen stations were used as part of the cross-validation (independent data set). Accordingly, the average density of the meteorological stations used to interpolate the grid-based map of rainfall erosivity (interpolation dataset, n = 370) was one in every 27.8 km × 27.8 km (or 775 km 2 ). The elevation of the stations varies from −3 to 2541 m a.s.l. (average 468 m a.s.l.), with a distribution following the altitudinal gradient of the country. Forty-two per cent of the total monitored stations are located in lowlands (0-300 m a.s.l.), whereas 34 and 24% are located in hilly (300-700 m a.s.l.) and mountainous areas (>700 m a.s.l.). The stations above 1000 m a.s.l. (12% of the total monitored stations) are primarily located in the Alpine region and the Northern Apennines.
The database consisted of 3408 rainfall years, with a mean of 8.8 years per station ( Table 1). The vast majority of the stations covered the same number of years of records, minimizing measurement errors related to the inter-annual variability of rainfall erosivity (Capolongo et al. 2008;Borrelli et al. 2014). The pluviographic series showed a low number of missing data (x ̅ = 0.4%), which were processed with a statistical procedure to ensure the high quality of the R-factor measurements.

Calculation of rainfall erosivity
Monthly rainfall-runoff erosivity values (R-factor) (MJ mm h −1 ha −1 yr −1 ) were computed following the RUSLE handbook instructions (Renard et al. 1997). According to Wischmeier and Smith (1978), the rainfall erosivity index (EI 30 ) of each erosive storm between January 2002 and December 2011 is given by the product of the storm kinetic energy and the maximum 30-min intensity (I 30 ): where e r is the unit rainfall energy (MJ ha −1 mm −1 ) and v r is the rainfall volume during a given period (r). The unit rainfall energy is calculated for each time interval according to Brown and Foster (1987): where I r is the rainfall intensity during the time interval expressed in mm h −1 . Finally, the annual average rainfall erosivity is expressed as where the (EI 30 ) i is the EI 30 for rainstorm i and j is the number of rainstorms in an N-year period. Storm-by-storm summaries of precipitation, duration, intensity, kinetic energy and the rainfall erosivity index (EI 30 ) were computed using the Rainfall Intensity Summarization Tool (RIST v3.88) (USDA 2014). A threshold of 1.27 mm of precipitation in six hours was selected to represent breaks in rainstorms. Afterwards, rainstorms with less than 12.7 mm of precipitation were omitted from the EI 30 calculation (Foster et al. 1981;Renard et al. 1997). An exception was made for rainstorms with at least one peak that is greater than 6.35 mm during a period of 15 min (Foster et al. 1981).
The vast majority of the pluviographic data were sub-hourly records (30-min time-step), comprising 91.6% of the entries. The remaining 8.4% of the entries, with a 60-min time-step, were located in Sicily, Ligura and Molise. According to the literature (Williams and Sheridan 1991;Yin et al. 2007), the R-factor (Wischmeier and Smith 1978) is underestimated as the time step increases from 5 to 60min. A linear interpolation was employed to homogenize the R-factor values obtained from two different sets of time resolution data (Yin et al. 2007;Capolongo et al. 2008;Panagos et al. 2015). The statistical relationships were assessed on a monthly basis by computing the EI 30 rainstorm values of neighbour stations using both 30-and 60-min time-step data.
With respect to the meteorological stations with 60-min time-step pluviographic records, the RIST (USDA 2014) outcomes for 37 (9.6%) stations were adjusted to avoid underestimations of the EI 30 and R-factor related to the different temporal resolutions. Linear regression analyses between 30-min and 60-min time-step data rainstorm erosivity showed monthly correction factors ranging from 1.15 to 1.41 (x ̅ = 1.27) in the North and from 1.19 to 1.47 (x ̅ = 1.32) in the South of Italy ( Table 2). The results showed a high coefficient of determination (R 2 averaging 0.972), consistent with previous findings (Yin et al. 2007;Panagos et al. 2015).
Inter-annual data gaps of the pluviographic records did not exceed 5%. However, monthly values of the R-factor (R m ) were statistically filled (Rudolf and Schneider 2005) to avoid even a small underestimation: where R mc is the computed monthly R-factor and MDC is the monthly data coverage ≥95%.

Spatial prediction of the R-factor
Given the high observation density and the rather similar intra-annual variability of rain patterns across the study area, a regression-kriging model (Odeh, McBratney, and Chittleborough 1995) based on a stepwise approach (Hengl, Heuvelink, and Stein 2004) was employed to interpolate the at-site estimated R-factor to a map. This approach made predictions by modelling the statistical relationship between the target variable (i.e. rainfall erosivity) and auxiliary environmental variables (covariates) (Goovaerts 1999) at sample locations. Subsequently, the target variable in the unsampled locations was predicted using the known value of the auxiliary variables (Hengl, Heuvelink, and Rossiter 2007). The model first employs regressions of auxiliary information and then employs simple kriging to interpolate the residuals from the regression model: wherem(s0)represents the deterministic part fitted by the regression model and1(s0)represents the interpolated residuals.
The regression-kriging model was designed and run using R language. The interpolation dataset was composed of average annual R-factor values of 370 meteorological stations. A set of 23 environmental covariates (e.g. climate data, land surface and spatial parameters) was created in an ArcGIS 10.2 environment and imported in the regression-kriging model. Covariates were tested for their significance using the stepwise regression method to identify the best (Meusburger et al. 2012) subset of predictors. For the predictors, a threshold for a significance of 0.1 was set. Accordingly, seven covariates were identified as being significant: (1) elevation, (2) annual rainfall, (3) average temperature of the coldest quarter, (4) isothermality, (5) temperature seasonality, (6) maximum temperature of the warmest month and (7) precipitation of the coldest quarter. A stable semivariogram model was employed to parameterize the spatial autocorrelation between the residuals derived from the multiple regression analysis.
The regression-kriging model performance was tested for both fitting and cross-validation. The amount of variation explained by the model was tested by leave-one-out cross validation (Cawley and Talbot 2004). A procedure consisting of fitting the model n-1 timeswith n being the number of observations in the datasetand at each time, one observation is left out of the fitting sample. A further assessment of the regression-kriging prediction accuracy was undertaken using an independent validation dataset (n = 48) composed of 16 stations that were randomly selected from our overall dataset (n = 386) and 32 stations that were derived from the literature (Diodato and Bellocchi 2010).

Comparison with approaches based on the rainfall volume
The results obtained in this study were compared with previous grid-based maps of rainfall erosivity proposed by Montanarella (1999,2000) (1 × 1 km cell size) and Bosco et al. (2015) (ca. 25 × 25 km cell size) for modelling soil erosion loss in Italy and Europe. Montanarella (1999, 2000) proposed a model where the average monthly rainfall is the only input required to estimate rainfall erosivity. A simple algorithm assumes that the R-factor is proportional to the cumulative yearly rainfall (Diodato and Bellocchi 2007). Bosco et al. (2015) combined auxiliary variables, covariates and seven empirical equations with the aim of improving the studies of Montanarella (1999, 2000). For the Italian region, the relationship obtained by Ferro, Porto, and Yu (1999) through a comprehensive rainfall erosivity analysis in the South of Italy was employed.
The purpose of this comparison was (1) to evaluate the performances of these previous methods based on rainfall volume and (2) to assess the spatial differences with respect to the approach proposed in this study.

Rainfall erosivity in Italy
In this study, 3408 years of continuous pluviographic records were processed using the RIST software (USDA 2014), identifying 70,095 erosive rainstorms. The average annual R-factor of the investigated meteorological stations totals 1715.2 MJ mm h −1 ha −1 yr −1 , with a standard deviation of 1171.2 MJ mm h −1 ha −1 yr −1 . This value is more than two times higher than the average R-value obtained by the 1290 non-Italian stations (723 MJ mm h −1 ha −1 yr −1 ) contained in the Rainfall Erosivity Database on the European Scale (REDES) . The average rainstorm erosivity index (EI 30 ) is 81.2 MJ mm h −1 ha −1 , with an average rainstorm duration of 20.9 h and an average maximum 30-min intensity of 12.2 mm h −1 .
The spatial distribution of the annual average rainfall erosivity (Figure 1) is highly variable within the country, with a 97% difference between the lowest values (216.9 MJ mm h −1 ha −1 yr −1 ) (Val di Breguzzo station, Trento) and the highest values (6939 MJ mm h −1 ha −1 yr −1 ) (Subit di Attimis, Friuli-Venezia Giulia). The north-westernmost region of the country (Friuli-Venezia Giulia) shows the highest average annual R-factor values. Here, half of the stations reported rainfall erosivity values above 4000 MJ mm h −1 ha −1 yr −1 , with an average of 4581 MJ mm h −1 ha −1 yr −1 . A situation consistent with the findings reported in the literature by Petan et al. (2010) was found in the neighbouring stations of western Slovenia in which the maximum European R-factor values were calculated . Subsequently, the areas with a higher frequency of thunderstorms (van Delden 2001), for example, Liguria, the northern border between Lombardy and Piedmont, Campania, Veneto, and to a lesser extent, Calabria and northern Tuscany, locally show high (2500-3500 MJ mm h −1 ha −1 yr −1 ) to very high (3500-5000 MJ mm h −1 ha −1 yr −1 ) rainfall erosivity. R-factor values below the national average were observed in the eastern sectors of the Apennine chain, with much of the Po Valley, Trentino, Valle d'Aosta, Basilicata and the major islands. Overall, 14% of the meteorological stations (n = 56) reported in this study show R-factor values below the European average. These stations are mostly concentrated in the driest sectors of the Alpine regions (Valle d'Aosta, Trentino-Alto Adige) and in Apennine intermountain valleys.
A cross-regional comparison of the descriptive statistics of the R-factor (NUTS-2 level; Eurostat 2014) is reported in Table 3. The rainfall erosivity was further evaluated in the context of climatic zones (EEA 2014) ( Table 3). The Mediterranean zone shows an average rainfall erosivity slightly higher than the other two main biogeographical zones (Alpine and Continental) (Figure 2). Seasonal data analysis did not reveal substantial differences in the R-factor distribution throughout the year between the Alpine and Continental biogeographical zones. In contrast, for the Mediterranean zone, a greater seasonality of the R-factor was observed, showing a more pronounced erosivity during the coldest months.
The average annual values of the erosive rainstorm ranged from 143 to 2677 mm per year. The monthly rainstorm amount was heavy during October, November and December, with maximum values occurring in November (Figure 3). The pluviometric regime decreased during the spring months, reaching minimum values in July. The number of erosive rainstorms per yearcomputed for the stations with 10-year record lengthsshows a bimodal distribution that was consistent with the general pluviometric regime. An average of 20 erosive rainstorms per year was observed, ranging from 6.7 at the station of Manfredonia (Puglia) to 35.1 at the station of Farra di Soligo (Veneto).
The temporal variability of the rainfall erosivity throughout the year shows high values during the autumn period and during summer. The minimum and maximum monthly values of the R-factor do not coincide temporally with the driest and wettest months (Figure 3). The highest average monthly rainfall erosivity is recorded in September (257.5 MJ mm h −1 ha −1 yr −1 ), whereas the lowest is recorded during February (48.9 MJ mm h −1 ha −1 yr −1 ). The rainstorms occurring during the summer months (July to September) (21.1% of the annual total rainstorms) account for more than onethird of the total R-factor (37.4%). During this season, the rainstorms are characterized by high intensity events (x ̅ maximum 30-min rainfall intensity 19.5 mm h −1 ) of short duration (x ̅ = 12.3 h).
The analysis of the rainstorm erosivity density (ED) (Kinnell 2010) for the 70,095 erosive events shows a high temporal variability (Figure 4), with values ranging from 0.88 MJ ha −1 h −1 (February) Table 3. R-factor descriptive statistics per region (NUTS-2) and biogeographical zone (*) (EEA 2014 to 4.36 MJ ha −1 h −1 (July). The erosivity density of the summer months is, on average, 2.5 times larger than for the rest of the year. Compared to the general European conditions ) (x ̅ = 0.92 MJ ha −1 h −1 ), Italian rainstorms are characterized by a considerably higher erosivity density (x ̅ = 2.01 MJ ha −1 h −1 ). The block-diagrams reported in Figure 5 display the spatio-temporal trends of the average seasonal R-factor values, offering a three-dimensional view of the spatio-temporal patterns described in the previous paragraphs. During the winter months (January to March), where rainfall erosivity is rather low, two trends of rainfall erosivity oriented West-East and North-South are observed ( Figure  5(a)). In this season, the southernmost regions (Calabria, Eastern Sicily and Apulia) show rainfall erosivity values higher than those of the rest of the country. However, this situation reverses moving towards the summer months, where a clear South-North-oriented trend occurs ( Figure 5(b)-(c)). An intermediate situation is noticeable during autumn, with a less remarkable North-South trend ( Figure 5(d)).
The analysis of the rainstorm EI 30 index showed that its spatio-temporal distribution was slightly different from the R-factor, showing high (80-100 MJ mm h −1 ha −1 ) to severe values (100-300 MJ mm h −1 ha −1 ) in areas where the R-factor was not very strong (800-1300 MJ mm h −1 ha −1 yr −1 ), that is, primarily located in the Mediterranean zone along the Ionian coasts, the southern Adriatic coasts and the two major islands. Conversely, some areas along the central Apennines, characterized by a higher R-factor (>1500 MJ mm h −1 ha −1 yr −1 ), showed considerably lower EI 30 values (55  mm h −1 ha −1 yr −1 ). The analysis of the pluviographic records for these locations, showing inverse trends of the EI 30 index and R-factor, revealed rather different rainstorm dynamics. In the first case, the rainstorms of the Mediterranean climate are characterized by high temporal variability and flashy characteristics (Angulo-Martínez et al. 2009;Fiori et al. 2014). Behaviour strongly affects the rainfall erosivity (González-Hidalgo, Peña-Monné, and de Luis 2007), which here depends on the occurrence of high-intensity rainstorms of short duration. In the second case, the lower erosivity of rainstorms is offset by the large number of erosive events occurring throughout the year. In the Apennine area, events mostly related to northwestern flows (Fiori et al. 2014) trigger orographic rainfall in winter and convective storms in summer (Miglietta and Rotunno 2010). Seasonal analysis of the EI 30 index shows two clear trends oriented West-East and North-South during the winter months (January to March) (Figure 6(a)). Similar but less pronounced trends are notable for the secondary cold season (October to December) (Figure 6(d)). During spring and summer (Figure 6(b-c)), the North-South trend previously observed is missing and the rainstorms present higher EI 30 index values.

Spatial modelling
In this study, the spatial distribution of the average annual rainfall-runoff erosivity (R-factor) over the Italian Peninsula was modelled through a regression-kriging approach. The result of the interpolation procedure is illustrated in Figure 7. The average annual rainfall erosivity totals 1351 MJ mm h −1 ha −1 yr −1 , with a range of 175 to 6985 MJ mm h −1 ha −1 yr −1 . As the definition of high rainfall erosivity depends on the study location , we adopted a statistical approach to classify the rainfall erosivity into 10 classes determined through a quantile-based analysis. Accordingly, 60% of the study area is classified by low to medium annual erosivity values (classes 1-6, 175-1430 MJ mm h −1 ha −1 yr −1 ). For the rest of the area (40%), the annual erosivity is classified as medium-strong (classes 7-8, 1430-1830 MJ mm h −1 ha −1 yr −1 ) and strong (classes 9-10, 1830-6985 MJ mm h −1 ha −1 yr −1 ).
The equation obtained using the multivariate linear model showed a good performance for both the leave-one-out cross-validation (R 2 cv = 0.777) and the fitting dataset (R 2 = 0.779). The annual R-factor was significantly related to elevation, annual rainfall, average temperature of the coldest quarter, isothermality, temperature seasonality, maximum temperature of the warmest month (α < 0.01), and to a lesser extent, precipitation of the coldest quarter (α < 0.05). The residuals derived from the multiple regression analysis were interpolated by an ordinary kriging using a semivariogram model. The predictions of the regression-kriging for the independent validation dataset (n = 48) showed an R 2 of 0.39 (Figure 8(a)). However, the scatterplot between the calculated and predicted R-factor values considering only the stations derived from this study (n = 16) illustrates a better correlation, with a high coefficient of determination (R 2 = 0.82) (Figure 8(b)).
The results presented in this study confirmed the complex spatio-temporal patterns of the R-factor observed in regional studies in Italy (Aronica and Ferro 1997;D'Asaro, D'Agostino, and Bagarello 2007;Capolongo et al. 2008). Despite this, the geostatistical approach proposed in this study to interpolate the R-factor provided satisfactory results that were in line with other studies conducted in Europe (Goovaerts 1999;Lloyd 2005;Hengl, Heuvelink, and Rossiter 2007

Comparison with approaches based on rainfall volume
The results of this study were compared with the maps of rainfall erosivity proposed in previous studies conducted by the Soil Team of the Joint Research Centre. The comparison results showed a weak correlation between the at-site average annual rainfall erosivity values measured in this study and those predicted by Montanarella (1999, 2000) (R 2 = 0.014) and Bosco et al. (2015) (R 2 = 0.199). Both maps tend to underestimate the annual R-factor values.   . Comparison between the rainfall erosivity map of Italy obtained in this study through a regression-kriging approach (500 × 500 m cell size) (left) and the maps proposed by Montanarella (1999, 2000) (1000 × 1000 m cell size) (middle) and by Bosco et al. (2015) (ca. 25 × 25 km cell size) (right). Rainfall erosivity values (R-factor) are ranked into five classes determined through quantile-based analyses. Montanarella (1999, 2000) showed an average underestimation of −53%, whereas for Bosco et al. (2015), the underestimation was −58%. This obviously leads to a low agreement between the maps proposed by Montanarella (1999, 2000) and Bosco et al. (2015) and that presented in this study. Although the three maps may share some spatial patterns of rainfall erosivity (Figure 9), the predicted values present marked differences in magnitude over the entire country area. The maps of Montanarella (1999, 2000) and Bosco et al. (2015) show average rainfall erosivity values of 809.6 (σ = 72.3) and 665.6 MJ mm h −1 ha −1 yr −1 (σ = 307.4), respectively. Therefore, these previous estimates were considerably lower than that achieved by the regression-kriging model (1351.2 MJ mm h −1 ha −1 yr −1 ; σ = 600.8).
The highest underestimation occurred in coincidence with the stations showing high rainfall erosivity, particularly those located in Veneto, Friuli-Venezia Giulia, Piedmont and Lombardy. Underestimations are also pronounced in southern Calabria, on the border between Lazio and Tuscany and eastern Sicily. Slight overestimations are observable in areas characterized by low rainfall erosivity of the Northern Apennines, Trentino-Alto Adige, Campania-Lucania Apennines and Central inland Sicily. From this analysis, it is inferable that the approximation equations based on the rainfall volume proposed by Montanarella (1999, 2000) and Bosco et al. (2015) cannot satisfactorily explain both the magnitude and patterns of rainfall erosivity in Italy.

Data availability
The rainfall station dataset presented in this study and part of the REDES project (Rainfall Erosivity Database on the European Scale; Panagos et al. 2015) is available for download on the European Soil Data Centre (ESDAC) web platform. The dataset is provided in an ESRI shapefile format and may encourage further studies into the spatial variability of rainfall erosivity in Italy.

Conclusions
The Universal Soil Loss Equation (USLE) and its revised versions (RUSLE, MUSLE, WaTEM/ SEDEM) are still the most used soil erosion prediction models in Italy. Therefore, a reliable estimation of the rainfall erosivity index is a crucial precondition for the successful identification of soil erosion risk areas in which soil conservation practices are needed. This analysis of rainfall erosivity yielded new information for soil erosion assessments in Italy. Using pluviographic records from 386 meteorological stations distributed across the country, new insights into the geography of rainfall erosivity in Italy were achieved, disclosing 'when' and 'where' rainfall hit the ground with higher erosive power.
A regression-kriging approach was successfully employed to map the rainfall erosivity in Italy at a 500-m spatial resolution. This achievement was completed by analysing 3408 years of continuous pluviographic records and more than 70,000 erosive rainstorms. The spatial interpolation model showed a very good performance (R 2 cv = 0.777 for the cross-validation, R 2 = 0.779 for the fitting). The annual rainfall erosivity values showed high spatial variability ranging from 170 to 6978 MJ mm h −1 ha −1 yr −1 . A large part of the study area (20%) revealed annual erosivity values greater than 1800 MJ mm h −1 ha −1 yr −1 . It was also shown that a considerable underestimation of rainfall erosivity characterizes the previous estimates reported in the literature ( Montanarella 1999, 2000;Grimm et al. 2003;Bosco et al. 2015). The new map presented in this study constitutes an important source of information for predicting soil erosion in Italy. Moreover, it has paved the path towards further investigations and interpolations of the seasonal spatio-temporal variability of rainfall erosivity.