Historical predictability of rainfall erosivity: a reconstruction for monitoring extremes over Northern Italy (1500–2019)

Erosive storms constitute a major natural hazard. They are frequently a source of erosional processes impacting the natural landscape with considerable economic consequences. Understanding the aggressiveness of storms (or rainfall erosivity) is essential for the awareness of environmental hazards as well as for knowledge of how to potentially control them. Reconstructing historical changes in rainfall erosivity is challenging as it requires continuous time-series of short-term rainfall events. Here, we present the first homogeneous environmental (1500–2019 CE) record, with the annual resolution, of storm aggressiveness for the Po River region, northern Italy, which is to date also the longest such time-series of erosivity in the world. To generate the annual erosivity time-series, we developed a model consistent with a sample (for 1981–2015 CE) of detailed Revised Universal Soil Loss Erosion-based data obtained for the study region. The modelled data show a noticeable descending trend in rainfall erosivity together with a limited inter-annual variability until ~1708, followed by a slowly increasing erosivity trend. This trend has continued until the present day, along with a larger inter-annual variability, likely associated with an increased occurrence of short-term, cyclone-related, extreme rainfall events. These findings call for the need of strengthening the environmental support capacity of the Po River landscape and beyond in the face of predicted future changing erosive storm patterns.


INTRODUCTION
Hydrological extremes [1][2][3] are an important component of climate variability and are partly governed by ocean dynamics 4 . Holocene landscapes are generally believed to be more or less in an equilibrium mode with external environmental driving forces [5][6][7] . However, moderate changes in storminess are accompanied by changes in hydrological extremes [8][9][10] . Storm (rainfall) erosivity, i.e. the capacity of rainfall to cause soil erosion, depends on these extremes. Greater awareness on the temporal variability of storm erosivity is important as it has implications for the understanding of geomorphological dynamics such as erosional soil degradation 11 and other landscapes stresses like flash-floods and surface landslides 12 . In fact, most of the damaging hydrological events worldwide are associated with strong erosive rainfall [13][14][15][16][17] . Erosive storm events usually result in accelerated soil erosion 18,19 and mud floods in urban areas 20,21 , and the consequences can be severe in terms of reduced socio-economic sustainability [22][23][24] .
An intensification of the hydrological cycle [25][26][27] , in particular, is likely to increase the magnitude of erosive events 28,29 as well as storm runoff extremes 30 . When the time-lapse decreases between successive shocks 31 , the responses of landscapes to changing disturbance regimes are also increasingly likely to be affected by the linkage of past damaging hydrological events 32 . In this way, shifts in erosive rainfall patterns could have a large impact due to a more powerful hydrological cycle and a concentration of storms in sporadic and more erratic and exacerbation events [33][34][35][36] . Despite uncertainties in total precipitation changes, extreme daily rainfall averaged over both dry and wet climate regimes show robust increases 15 in parts of the world. In particular, the time-variability of storminess and its unequal distribution create potential changes in the power of rainfall to drive erosional processes in a landscape 37,38 .
Advances have been made in the understanding of the dynamics of past and future extreme precipitation worldwide 26,39 . However, while atmospheric thermodynamics may suffice to explain how precipitation is started, the deterministic forecast of the intensity and spatial extent of storms is often inaccurate. This is mainly a result of the intrinsic uncertainty of precipitation, which should also be based on probabilistic, statistical and stochastic descriptions 40,41 . For the past, the reconstruction of long timeseries of climate extremes supports that naturally climatic processes are governing the occurrence of erosive rainfall and, thus, help us to understand present processes and to improve future projections 42 . Attempts have been made to relate weather types or climate indices, to temporally coincident multi-secular trends in hydrological processes 43,44 . In this way, the predictive performance of the trends leaves a small potential for long-term predictability in rainfall 45 . Recognising trends may be helpful in achieving mitigation and adaptation strategies to counter the possible consequences of abrupt changes in the frequency or severity of climate extremes. In such respect, research in historical climatology adds an additional dimension to state-of-the-art climate model simulation approaches of present-day and future predicted conditions. Simulations with state-of-the-art climate models are useful in studying climate changes, but they show considerable uncertainty in terms of internal hydroclimate variability 46,47 while exacerbating risks associated with extreme events 48 , especially for global-to-local estimation of erosion 49,50 . An approach to address the uncertainty produced by climate models is through the so-called stochastic synthesis, based on the function of auto-covariance of different hydroclimatic processes 51 , which directly stimulates the variability, stochasticity and fractality that seems to be inherent and consistent in geophysical processes 41 .
Long documentary records can help relating local and regional hydrological extreme features to impacts 52,53 and climatic forcing 54,55 , and in integrated multi-disciplinary approaches to climatic variability in hazard-exposed areas 56,57 . Studies in historical climatology can contribute with valuable new perspectives to the research of landscape preservation in the light of climate change, especially in terrestrial river systems that are highly dynamic and regulated by intricate climatic, geomorphic and ecological processes 58 . The Po River landscape (PRL; Fig. 1), whose Po Delta Regional Park was designated in 1999 as a World Heritage Site by the UNESCO, is the area of Italy where the most damaging hydrological extreme events are recorded 59 , and where the population density of the country is the highest 60 .
The Po River Landscape (~70,700 km 2 ; Fig. 2A) extends over large portions of northern Italy, in the transition zone between Central Europe to the north and the Mediterranean region to the south. Small portions extend into Switzerland (5.2%) and France (0.2%). Geographic features change across the three main sectors of the landscape 61 . The Alpine sector (45% of the total area) forms an imposing mountain range to the north, west and south-west, with a length of~570 km and a width varying from 25 to 110 km. The Apennine sector (15% of the total area) consists of western hilly reliefs and minor reliefs bordering it to the south, with a length of~210 km and a width varying from 25 to 90 km. The Po River Valley (40% of the total area) develops as a strip of territory that extends along with the two previous sectors over~490 km and is 20 to 120 km in width. The basin has a roughly rectangular  Fig. 1 Perspective of the Po River landscape's hazard-exposed to the storms aggressiveness. Perspective of the Po River landscape's hazard-exposed to the storms aggressiveness, which develops as a strip of territory inserted between Alps (to the north and west) and Apennine (to the south); the winding curve (in dark blue) is the Po River that within the town of Pavia and surrounded cities and villages goes up towards the town of Turin and thus the Alps (arranged from OpenStreet Map https://demo.f4map.com/#camera. theta=0.9). 0 1000 500 Kilometers  123 .
shape and with a preferential east-west orientation of the surface drainage. The hydrographic network is set according to a dendritic drainage pattern. The highest point is the summit of Mont Blanc (4810 m a.s.l.), at the border between Italy and France, with an average altitude of~740 m. The climate of the study area is rather complex and diversified given its geographical position and the different morphology of the sectors composing it. The meteorological conditions of the PRL largely depend on the fronts formed between polar and tropical air masses. In particular, different situations are produced that depend on the cell that is interposed between the subpolar cold air and the warm Mediterranean air on one side, and between the damp maritime climate of the west and the dry continental climate of the east on the other 62 . In this way, the PRL climate shows a transitional character between the subcontinental climate of Central Europe (Alpine and Boreal) and the Mediterranean one (Warm Alpine) 63 . The average rainfall over the PRL is~1100 mm yr −1 , while for the mountain and Apennine sectors the precipitation is~2000 mm yr −1 , and for the lowland sector it is~750 mm (Fig. 2B). The precipitation regime is characterised by summer and early autumn maxima (mostly erosive), and winter and spring minima. The spatial pattern of mean annual of rainfall erosivity over the PRL shows considerable values in the central sector (~2000 MJ mm ha −1 h −1 yr −1 ), and until 4000 MJ mm ha −1 h −1 yr −1 , across Valle d'Aosta, northern Piedmont and the central Alps. In the western Alps, and in some limited areas over the low river basin, the erosivity is somewhat low, reaching values between 700 and 1000 MJ mm ha −1 h −1 yr −1 (Fig. 2C).
To obtain actual rainfall erosivity values following the Revised Universal Soil Loss Equation ((R)USLE) methodology 64 , careful precipitation measurements on short time-scales are required. In its original formulation 65 , the rainfall erosivity factor (R) is calculated by an empirical relationship of the measured 30-min rainfall intensities (R = E·I 30 ). For a given site, the monthly values of R (R m ) are given by the sum of all the single-storm EI 30 values for the month of interest. The term EI 30 (MJ mm h −1 ha −1 ) is the product of storm kinetic energy (E), calculated over time steps of minutes of constant storm intensity, and the maximum 30 minintensity (I 30 ). This makes it challenging to perform historical studies since measurements of this type are typically not available prior to the modern instrumental period. Hence, the only possibility is offered by modelling approaches, making use of long-term precipitation inputs. Model-based erosivity estimates back in time can consequently only be performed at a lower resolution, both in time (i.e. annual resolution or lower) and space (i.e. non-locally calibrated) 66 . For that, parsimonious models can be used since they overcome the limitations imposed by detailed models, which are data demanding and inapplicable to historical times 67 . In the event satisfactory instrumental input data are unavailable, data recorded in historical documentary sources can be used to support rainfall erosivity estimates 68 .
In this study, we first acquired comprehensive knowledge about potential drivers of rainfall erosivity in the PRL, where highintensity rainfalls occur from June to October due to the prevalence of thunderstorms during this part of the year in the region. Autumn precipitation is also triggered by synoptic disturbances coming from the North Atlantic Ocean, fed by moisture flukes from the hot surface of the Mediterranean Sea and maintained by mesoscale processes 69 . We use these factors (amounts of precipitation and weather anomalies) to develop a parsimonious model for reconstructing annual rainfall erosivity over the period 1500-2019 CE, allowing us to capture a wide range of climate variability, and identifying landscape-stress changes. The results of this study complement, for the past five centuries (1500-2019 CE), the twelve century-long (800-2018 CE) reconstructions of extreme hydro-meteorological events across Italy 56 and in the Po River Basin 70 .

Rainfall erosivity model calibration and validation
For the reconstruction of annually resolved storm erosivity (MJ mm ha −1 h −1 yr −1 ) for the PRL, we developed and calibrated the Rainfall Erosivity Historical Model (REHM, Eq. (1)), based on summer and autumn precipitation data and Gaussian-filtered values of the reconstructed annual severity storm index sum (ASSIS(GF)) from Diodato et al. 56 . For the calibration period 1993-2015 CE, we obtained the coefficients (2) and (4). With these values, the ANOVA returned a highly significant relationship (p~0.00) between observed and predicted erosivity values. The R 2 statistic (goodness of fit) indicates that the REHM explains 74% of the erosivity variability. p-value < 0.05 means that there is still a statistically significant relationship between the estimated and actual data. Though the R 2 statistic indicates that the model only explains 50% of the variability of actual erosivity, the time-variability as a whole is well reproduced (Fig. 3C). The MAE and the KGE are equal to 290 MJ mm ha −1 h −1 yr −1 and 0.66, respectively. Summary statistics (Table 1) for the overall period of available data (1981 show that the modelled data tend to underestimate the observed values, with higher means than medians indicating that in both modelled and actual time-series a few years had a much higher erosivity than most of the rest. This underestimation indicates that the model can only capture a part of the extreme values, thus lowering both the average and the respective percentiles. Whether the paired Student t-test (t = 2.09) for the difference between means indicates a marginal significance (p = 0.044), the Kolmogorov-Smirnov statistic (D = 0.23) for the difference between distributions did not show any significant difference (p = 0.320) between the modelled and observed time-series. The variability is also not statistically different between the observed and modelled time-series (variance ratio F-test = 1.27, p = 0.052).
Comparison to alternative models In order to evaluate how much simpler equations are able to estimate rainfall erosivity, we recalibrated three alternative models: a 2nd-order polynomial regression 67 : Z = 1178 − 1.263 · P ann + 0.00159 · P ann 2 (Fig. 4A), a logarithmic regression: Z = 1974 · LN(P sum + P aut ) − 10,877 (Fig. 4B), and a multiple linear regression model: Z = 2.886 · (1.400 · P sum + P aut ) − 319 (Fig. 4C), where Z represents the actual annual rainfall erosivity, P ann the annual precipitation, P sum and P aut , the summer and autumn precipitation, respectively.
The scatter-plots in Fig. 4 present a moderate fitting, with R 2 around 0.5-0.6, below the value of 0.74 obtained with the RHEM. Overall, we can conclude that simpler models are less adequate than RHEM to estimate rainfall erosivity. In particular, the insertion of ASSIS(GF) as predictor in Eq. (1) brings a significant improvement in the estimate of annual erosivity, which is difficult to estimate with simpler solutions.
Long-term reconstruction of rainfall erosivity in the PRL Though consistent with the general trend of actual erosivity data, which are representative of the fluctuations of rainfall erosivity in the PRL, model estimates may not always be adjusted to peaks of erosivity (e.g. high percentiles in Table 1). It can thus be assumed that also in past times high erosivity values may have occurred, which are not entirely captured in the long-term modelled timeseries. This may be related to non-linear dynamics implied by the Hurst (H) phenomenon 71 , i.e. the rate of chaos, with H > 0.5 corresponding to persistent Brownian motions (0.5 corresponding to a random process). The estimated H exponent (R/S method 72 ) equal to 0.84 (above the threshold of 0.65 used by Quian and Rasheed 73 to identify a series that can be predicted accurately) reflects the existence of low-and high-intensity storm clusters in the erosivity series, due to the non-linear climate dynamics governing the occurrence of rainfall extremes 74 . However, the similar distribution and variability to those of the actual data (as identified during the observational period) support the use of model estimates to study trends and fluctuations even though in reality the erosivity peaks may be higher than the modelled ones. Figure 5 shows the rainfall erosivity evolution over the period 1500-2019 CE. After this reconstruction by means of Eq. (1), the time-series was analysed to find out possible patterns of rainfall aggressiveness and to compare contemporary with historical rainfall erosivity anomalies and their variability. The mean values of rainfall erosivity over the entire series is 1844 CE MJ mm ha −1 h −1 yr −1 . Anyhow, the temporal evolution of annual values presents a decreasing trend (Fig. 5A, blue curve and related polynomial), with change points detected in the years 1657 CE (statistics of Buishand 75 , and Pettitt 76 ) and 1709 CE (statistics of Alexandersson 77 , and Worsley 78 ), during the grand solar (Maunder) minimum of 1645-1715 CE 79 . We refer to 1709 CE hereafter, as the starting point of a calmer but more changeable period (Fig. 5A). However, stormy years (with greater erosivity values than the 98th percentile) may occur during any time period.
In this way, the landscape was subject to considerable stress either before or at the start of the Maunder minimum (e.g. 1529, 1579 and 1647 CE) and subsequently (1810, 1976 and 2019 CE; this last with blue dot in Fig. 5A), suggesting that the pressure exerted by aggressive rainfalls can return at periodical intervals. The wavelet power spectrum mostly indicates a~22-year periodicity (Fig. 5B), which reflects the~22-year magnetic polarity cycle of sunspot activity 80 . Similarly, Zanchettin et al. 81 found solar-type periodicities, suggesting that the sun may be one of the precursors of hydrological processes in northern Italy, involving the magnetic activity of sunspots.  The periods before and after 1709 CE reveals no similar patterns. Average values indicate that the period 1500-1708 CE, with 2098 MJ mm ha −1 h −1 yr −1 , was characterised by more erosive events. The second period, 1709-2019 CE, was instead affected by a less marked regime of erosivity, with 1682 MJ mm ha −1 h −1 yr −1 on average, but accompanied by a more changeable variability pattern (orange curve in Fig. 5A). While a quasiparallel trend of erosivity (blue bold curve) and coefficient of variation values is visible until 1708 CE, the two curves tend to divergence after this date, the coefficient of variation indicating a kind of exponential growth, passing from 0.20 to 0.35 between 1709 and 2019 CE (orange curve in Fig. 5A).
This means that the PRL has been exposed to progressively increasing inter-annual variability of rainfall erosivity, which translates into greater hazardous landscape stress. As a consequence, disasteraffected areas have become more vulnerable to climate because rainfall aggressiveness has become more changeable. This variability may be caused by different, time-scale-dependent processes (e.g. von Storch and Zwiers 82 ), with erosivity records showing rapid step-like shifts in variability occurring over decadal and multi-decadal timescales. In particular, extremes may persist and evolve in a more unexpected way. In fact, during recent decades, the aggressiveness seems to indicate a recovery in the erosive activity of rains. This trend is in accordance with a marked increase of erosivity density (the ratio of rainfall erosivity to precipitation, data not shown), which confirms an increased inter-annual variability. Though the frequency of occurrence of daily precipitation was found to decrease over in Italy 83 , episodes with shorter duration (from 1 to 3 h) have instead enhanced the torrential character of seasonal rains 83,84 . This accords with what reported by Colarieti Tosti 85 , which provided sufficient data to show that in the coming decades the polar vortex will undergo a phase of expansion towards the south with consequent exacerbation of the hydrological cycle in the central-western Mediterranean area. Already in 2019, torrential downpours and extreme flooding have battered many parts of Northern Italy, e.g. Venice (northeast) saw record-breaking water levels and the worst flooding in 50 years, and a viaduct near Savona (northwest) was washed away by a mudslide (https://tinyurl.com/tt7apoq).
In northern Italy, short-duration extreme hydrological events have risen over the last 70 years 86,87 . During 2004-2016, the highest number of such extreme events was mostly recorded in the central Piedmont region 88 . This trend is in agreement with regional climate simulations predicting increasing precipitation over the high Alpine elevations in the coming decades, associated with increased convective rainfall due to enhanced potential instability by high-elevation surface heating and moistening 89 . An increase in the frequency of damaging convective weather events over hazard-exposed landscapes is also expected over much of Europe towards the end of the twenty-first century 90 .
Climate model simulations with the new generation of Coordinated Downscaling Experiment over Europe (EURO-COR-DEX 91 ) reveal an increasing trend towards a higher frequency of river floods extremes across most of Europe with future global warming. Towards the end of the twenty-first century, the results from the EURO-CORDEX simulations show no significant change in annual precipitation over northern Italy, but at the same time an approximately 20% increase in maximum daily precipitation 92 . In order to assess whether there was a change in the frequency distribution, before and after the change-point of 1709 CE, we calculated the frequency histograms pattern for both these periods. It results in that before the change-point (Fig. 6A), the distribution of the erosivity was approximately normal, while later (Fig. 6B), the pattern appears to have a quasi-Poisson distribution.
The mesoscale atmospheric circulation can evolve into landscapeimpact weather systems such as rainfall aggressiveness and related damaging hydrological events. We improved the understanding of hydrological extremes, and in turn rainfall erosivity variability, thanks to the long continuous time-series (1500-2019 CE) developed for the PRL from climate records. This was performed with an integrated approach to reconstruct and assess rainfall erosivity for individual years. It included the collection and evaluation of climate documentary sources and a rainfall erosivity model for the reconstruction of past climate erosive variability and detection of climate signals. For the multi-centennial period assessed in this work, the main conclusions can be summarised as follows: 1. Annual rainfall erosivity across the PRL from 1500 CE onwards shows alternate stormy and quiet periods, which may be partially explained by a long-range dependence and non-linear climate dynamics of the pulsed occurrence of rainfall extremes. In this time-scale, a change-point is manifested during the Maunder minimum of solar activity (1645-1715 CE), with remarkable differences before (1500-1708 CE stormy period) and after (1709-2019 CE with a calmer but more changeable period). 2. The generally cold period 1500-1708 CE appears characterised by very stormy years (1529, 1560, 1579, 1610, 1679 and 1703 CE). The following period, until 1809 CE, is a very quiet phase, in the transition to increased storminess activity while inter-annual variability is growing. 3. Since~1810 rainfall erosivity shows substantially dissimilar dynamics than in the previous periods: precipitation aggressiveness remains quasi-stationary, but some sparse storm erosivity events occur, such as those around the years 1953, 1960, 2001 and, 2014 CE. The way how in the past decades Mediterranean cyclones have been producing trends of rainfall extremes remains elusive 93 , with significant negative trends of cyclone frequency in spring often compensated by positive trends in summer 94 . Raible et al. 95 supports an increase of extreme cyclone-related precipitation, purely thermodynamically driven by the temperature increase and the Clausius-Clapeyron relation 96 . An intensification of geohydrological hazard owing to an increased occurrence of severe rain events during the past three to five decades is acknowledged for north-western Italy 97 and inner hilly areas of central Italy 98 . However, no sustained long-term trend was observed from 1760 CE over northern Italy (including the Northern Adriatic Sea) in spite of pronounced inter-annual and inter-decadal variability of storminess 99 .
It can be remarked that the present analysis at an annual scale may mask important variations manifesting at finer time-scales. The methodology applied in this study is, therefore, most appropriate to address inter-annual and inter-decadal time-scales, but does not capture daily to seasonal changes which may impact on hydrology and land conditions and different spatial scales. However, the presented continuous (provisional) rainfall erosivity reconstructions over the PRL are of great importance for all kind of climatic analyses dealing with natural variabilities at decadal and centennial time-scales. They can be used to study the low and high-frequency variability, and the characteristics and extremes of climate at sub-regional scale, and can be compared with modelproduced simulations of natural and forced (external and internal) variability for the past centuries. The results also imply that environmental management can use data from long historical time-series as a reference for decision making. For this period, an annual series of erosivity data for the PRL was obtained by averaging the annual erosivity values determined in each station located in the basin, and this dataset was used to determine model parameters (calibration dataset). Prior to 1993, only data from the western part of the PRL were provided by Piedmont Region stations since 1981 100 . Such (R)USLE data from the period 1981-1992 were used for an independent evaluation of model estimates (validation dataset). The following adjusted erosivity data-AdjRðUSLEÞ i -were obtained from the original R(USLE) data of the set of stations located to the west side of the Po River Basin, which is basically an extrapolation to the whole of the PRL:

METHODS
; 996, where the actual erosivity data were summed over the n locations (l) and linearly combined with average autumn (P Autl ) and summer (P Suml ) precipitation in each year (i). The coefficients of the regression Y = AdjRðUSLEÞ i versus X = LN P n l¼1 ðRÞUSLE i þ P Autl þ P Suml À Á were obtained for the period 1993-2015, for which basin-wide (R)USLE data were available (R 2 = 0.87).
(R)USLE model 103,104 appears in the following form: where n is the number of recorded years; m j is the number of erosive events during a given month j; k is the index of a kth single event; v r is the rainfall volume (mm) during the rth period of a storm, which splits into m parts; I 30 is the maximum 30-min rainfall intensity (mm h −1 ); i r is the rainfall intensity during the time interval (mm h −1 ); P j is the rainfall amount (mm) during a given month j.

Numerical and categorical dependent variables
The reconstruction of annually resolved storm erosivity for the PRL was based on the seasonal precipitation data for a grid covering most of Europe (30W-40°E/30°N-71°N) from 1500 to 2000 on 0.5°by 0.5°r esolution 105 . The CRU Global Climate Dataset 106 provides an extension of the seasonal precipitation dataset until 2019. A categorical variable, the annual severity storm index sum (ASSIS), was derived from Diodato et al. 56 to overcome the lack of historical information about rain intensity. ASSIS was derived from several sources by transforming documentary information into a record set to 0 (normal event), 1 (stormy event), 2 (very stormy event), 3 (great stormy event) and 4 (extraordinary stormy event).

Historical modelling of rainfall erosivity
For the historical reconstruction of annual rainfall erosivity (MJ mm ha −1 h −1 yr −1 ), we developed the Rainfall Erosivity Historical Model (REHM), which uses as inputs summer and autumn precipitation data (mm) and Gaussian-filtered values of the reconstructed ASSIS: where the terms φ and ψ are expressed in mm −1 while the scale parameter A converts the result between brackets in MJ mm ha −1 h −1 yr −1 . The
process term φ is: where α (mm −1 ) and k (mm −1 ) are empirical parameters. The term φ modulates the autumn erosive precipitation. It equals α when ASSIS(GF) = 0. Since floods occur generally during the autumn season, the Gaussian-filtered annual storm-severity index sum, ASSIS(GF), was selected as an input covariable to capture storm-energy during this season. The interaction between autumn precipitation and ASSIS(GF) interprets the fact that wet autumn driven by severe storms increases the risk of erosivity. The multiplicative component φ • P Aut supports the non-linear dependence of rainfall erosivity on precipitation intensity (e.g. D'Odorico et al.) 107 . The ASSIS(GF) values from the original ASSIS time-series were smoothed by applying the low-pass Gaussian filtering technique described by Gjelten et al. 108 . Weighting coefficients (w ij ) were applied to derive the Gaussian function, ASSIS(GF), for each year j: where x i is the data point at year i, σ is the standard deviation, and n is the number of years in the series. The series of filtered values is established by letting j run through all data points. To remove variations at shorter time-scales than approximately 10 data points in the time-series, 11-year window (ij) and σ = 3 years were chosen (e.g. Harris et al. 109 and Manara et al. 110 ), which is a way to modulate the impact on rainfall erosivity of the decadal variability and long-term trend of storm-severity. Summer precipitation is characterised by an intrinsic high variability. The intensity of summer rainfall is modulated by Eq.
(4), a modified version of the variation coefficient developed by Aronica and Ferro 111 : where SD is the inter-seasonal standard deviation and Max(Ps) is the maximum seasonal precipitation per year (mm). Winter and spring rainfalls have not been taken into account by the model because rainfall intensity is generally low in those seasons in northern Italy, with limited erosive force. The concept of the model is summarised in Fig. 7 following Diodato and Bellocchi 112 . It shows that in Northern Italy brief and intense convective storms are more common in summer while autumn rains (of long duration and low intensity) usually originate from broad mid-latitude frontal activity. They carry high volumes of rainwater thanks to orographically enhanced stratiform precipitation (expressed by high percentiles), which also cause large-scale hydrological processes like flooding 113 . These processes are captured in the Rainfall Erosivity Historical Model by the storm-severity index of Eq. (2). Highly intensive convective rainfall results instead in the splash-erosivity process interpreted by the variation coefficient of Eq. (4). In this case, raindrops with high kinetic energy may cause the local detachment of soil particles 114 . To assess the model, statistical analyses were performed with STATGRAPHICS 115 , with the graphical support of WESSA 116 and CurveExpert routines 117 . The mean absolute error (MAE, MJ mm ha −1 h −1 yr −1 ) was used to quantify the differences between actual and predicted erosivity values, while determination coefficient (0 ≤ R 2 ≤ 1, optimum) and slope (b = 1, optimum) assessed the linear relationship between the two series. The Kling-Gupta index (−∞ < KGE ≤ 1 118 ) was used as an efficiency measure, with KGE > −0.41 indicating that a model improves upon the means of observations as a benchmark predictor 119 . The Durbin-Watson statistic 120,121 was performed to test for auto-correlated residuals because large temporal dependence may induce spurious correlations 122 . ANOVA p-values were used to present the statistical significance of the regression between estimates and the actual data. Student t-test and Kolmogorov-Smirnov test were applied to determine whether the two groups (estimates and actual data) originate from the same population. The variance ratio F-test was used to test the difference between variances.