==== Front Sci Rep Sci Rep Scientific Reports 2045-2322 Nature Publishing Group UK London 78857 10.1038/s41598-020-78857-3 Article Fingerprint of climate change in precipitation aggressiveness across the central Mediterranean (Italian) area Diodato Nazzareno 1 Ljungqvist Fredrik Charpentier fredrik.c.l@historia.su.se 234 Bellocchi Gianni 15 1 Met European Research Observatory, International Affiliates Program of the University Corporation for Atmospheric Research, Via Monte Pino snc, 82100 Benevento, Italy 2 grid.10548.380000 0004 1936 9377Department of History, Stockholm University, 106 91, Stockholm, Sweden 3 grid.10548.380000 0004 1936 9377Bolin Centre for Climate Research, Stockholm University, 106 91, Stockholm, Sweden 4 grid.462826.c0000 0004 5373 8869Swedish Collegium for Advanced Study, Linneanum, Thunbergsvägen 2, 752 38 Uppsala, Sweden 5 grid.494717.80000000115480420Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, 63000 Clermont-Ferrand, France 16 12 2020 16 12 2020 2020 10 2206216 6 2020 1 12 2020 © The Author(s) 2020Open 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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by/4.0/.Rainfall erosivity and its derivative, erosivity density (ED, i.e., the erosivity per unit of rain), is a main driver of considerable environmental damages and economic losses worldwide. This study is the first to investigate the interannual variability, and return periods, of both rainfall erosivity and ED over the Mediterranean for the period 1680–2019. By capturing the relationship between seasonal rainfall, its variability, and recorded hydrological extremes in documentary data consistent with a sample (1981–2015) of detailed Revised Universal Soil Loss Erosion-based data, we show a noticeable decreasing trend of rainfall erosivity since about 1838. However, the 30-year return period of ED values indicates a positive long-term trend, in tandem with the resurgence of very wet days (> 95th percentile) and the erosive activity of rains during the past two decades. A possible fingerprint of recent warming is the occurrence of prolonged wet spells in apparently more erratic and unexpected ways. Subject terms Climate sciencesEnvironmental sciencesHydrologyhttp://dx.doi.org/10.13039/501100004359Vetenskapsrådet2018-01272Ljungqvist Fredrik Charpentier Stockholm UniversityOpen Access funding provided by Stockholm University issue-copyright-statement© The Author(s) 2020 ==== Body Introduction Both natural and anthropogenic climate change can alter storm (rainfall) erosivity or the R-factor, i.e., the power of rainfall to cause soil erosion as defined in the Universal Soil Loss equation1 and updated versions of it2–4, due to the modification of rainfall patterns in time and space5. The oscillatory behaviour of extreme rainfall variability, and its associate erosive power, have huge impacts on agriculture6, hydrological processes7 and socio-economic dynamics8 across multiple spatial and temporal scales9, 10. Advances have been made in recent scholarship towards the understanding of the dynamics of past and future extreme precipitation worldwide11, 12. However, traditional climate extreme indices and large-scale multi-model inter-comparison studies, used for future projections of extreme events and associated impacts, often fall short in capturing the full complexity of impact systems13. Simulations with state-of-the-art climate models show noticeable uncertainty in terms of internal climate variability14 and climate response15, 16. In addition, only low temporal resolution (e.g., monthly and seasonal) precipitation time-series are available for longer time-periods, which reflect our still incomplete knowledge, and an inability to reconstruct long time-series back in time of the spatial and temporal distribution, as well as the magnitude and driving mechanisms, of the R-factor17. Rainfall erosivity is not only important for the understanding of surface-process dynamics such as erosional soil degradation17–19, and other landscape stressors like flash-floods and landslides20, but its dynamics also offer an opportunity to detect the fingerprint of recent climate change—especially in the Mediterranean region—which is a particular sensitive region regarding climate variability and a “hotspot” of climate change21. For parts of this region, the ongoing trend towards more extreme precipitation is expected to continue in the coming decades, contributing to the uncertainty in the projected occurrence and intensity of extremes22. High-resolution and well-dated records are needed to understand the long-term hydroclimatic variability in this region, but the careful rainfall measurements on sub-hourly time-scales, which are necessary to obtain actual rainfall erosivity values according to the (R)USLE (Revised Universal Soil Loss Equation) methodology4, are not available prior to the digital instrumental period starting in the 1980s23. In the absence of such detailed rainfall information, rainfall variability and its linkage to atmospheric pressure systems can be inferred from meteorological observations available for the Mediterranean region over the twentieth century24and even earlier back to mid-seventeenth century on a restricted regional basis25, but also from time-series of millennium-long hydrological extremes derived from documentary sources26. Long series derived from documentary data can help linking local and regional hydrological extremes to impacts27, assess climatic forcing features28, and recognise the climatic variability in hazard-exposed areas29, 30. Their analysis supports that climate can vary in response to natural processes such as those governing the occurrence of erosive rainfall and, thus, help us to understand present-day hydrological dynamics and improve projections of future changes31. In the Mediterranean region, local to regional climate experiences a mix of gradual and abrupt shifts of precipitating systems, including erratically distributed erosive storms32. In fact, while northern Europe is primarily influenced by the North Atlantic storm tracks33, southern Europe is situated at the junction of major air-mass circulations, periodically arising from Atlantic, Mediterranean and Siberia influences34. Different meteorological factors can contribute to the growth and recurrence of extreme rainfall, but two components appear essential for generating heavy precipitation over the Mediterranean region: a high water vapour content in the atmosphere and triggering events originating from thermodynamic or dynamic processes35. Deepening over the warm waters Mediterranean cyclones mostly form around a few centres, with a dominating region in the Gulf of Genoa, where a slowly moving low pressure field (or Vb-weather pattern) can bring large amounts of rainfall36. Modelling approaches making use of low-resolution precipitation data and documentary records of extreme weather events provide a means to derive long-term reconstructions. When satisfactory instrumental input data are unavailable, information from historical documentary sources can be used to support low-resolution storm-erosivity estimates37, both in time (i.e., with annual resolution or finer) and space (i.e., non-locally calibrated)38. In this study, we developed a modelling approach of the temporal fluctuations of storm-erosivity across the Mediterranean Central Area (MedCA), which is the most exposed region in southern Europe to aggressive rainfall (Fig. 1a,b). In this Mediterranean sub-region, high-intensity rainfall events can occur in any month due to the prevalence of thunderstorms throughout the whole of the year, which in turn can drive noticeable annual rainfall erosivity (Fig. 1c). On average for the period 2003–2012, rainfall erosivity over the MedCA ranged from over about 400 and 4000 MJ mm hm−2 h−1 year−1, with the western Tyrrhenian coast and the north-west being the most erosive-prone sectors with about 1500–2000 MJ mm hm−2 h−1 year−1 (Fig. 1c). In the eastern Alps and along the middle and low Adriatic (east) versant, the erosivity is somewhat lower, reaching values between 600 and 1000 MJ mm hm−2 h−1 year−1, respectively (Fig. 1c). Precipitation aggressiveness is generally triggered by synoptic disturbances coming from Atlantic cyclonic westerlies39, fed by heat and moisture fluxes from the Mediterranean Sea and maintained by mesoscale processes which, in turn, determine a suite of distinct, spatially-scattered erosive patterns (as in the schematic representation of Fig. 2).Figure 1 Study area. (a) Environmental setting. (b) Spatial pattern of mean annual rainfall erosivity (MARE) over central-southern Europe. (c) Map of MARE over Mediterranean Central Area for the period 1994–2013 (arranged with Geostatistical Analytics by ArcGIS-ESRI on the ESDAC-dataset, https://esdac.jrc.ec.europa.eu/content/global-rainfall-erosivity) 55. Figure 2 Scheme depicting the erosive rainfall pattern with cyclonic westerlies across a longitudinal transect west–east up to the middle Apennine of Italy. The city of Naples is shown, situated on the western coast of southern Italy (40° 50′ N, 14° 15′ E). Hillslopes are the dominant landform features in this area, where complex interactions between atmosphere and landscape systems operate at a range of time and spatial scales40. It follows that the mix of rain-producing erosivity events appears to vary not only spatially but temporally as well. However, detailed studies on the control factors of temporal dynamics are still limited, hence motivating the need of this study. In fact, some previous studies have determined rainfall erosivity using either limited rainfall records41–43 or non-uniform time intervals across different spatial scales44. Here, we used amounts of seasonal precipitation and weather anomalies to develop a simple, climatically interpretable model for reconstructing annual erosivity data in the MedCA over the period 1680–2019 CE. For both its length and spatial extension, this provided an unprecedented time-series to discern the climate change fingerprint and communicate the regional erosivity hazard in Central Mediterranean. In this study, we have excluded all data prior to 1680 because precipitation data for earlier times have been reconstructed from non-instrumental data exclusively. They are, thus, affected by a larger uncertainty than more recent instrumental data, especially for the southern sector of the MedCA45. Results and discussion Rainfall erosivity model calibration and validation To estimate areal-mean annual rainfall erosivity over the MedCA, a simplified statistical model was developed, which summarises the relationship between spatial patterns of climate and storm erosivity, consistent with a sample calibration (1994–2015 CE) and validation (1981–1993 CE) of detailed (R)USLE-based data obtained for the study area. We calibrated the Rainfall Erosivity Mediterranean Model (REMM, Eq. (1) in Methods), based on multiple observational time-scales, such as: seasonal precipitation, its variability and the Gaussian-filtered annual severity storm index sum SSI(GF) from Diodato et al.26. For the calibration period 1994–2015 CE, we obtained the coefficients A = 0.0512 MJ mm−1 hm−2 h−1 year−1, B = 0.1128 MJ hm−2 h−1 year−1, C = 637 MJ hm−2 h−1 year−1, α = 10.0 and k = 6.0 in Eq. (1). With these values, an ANOVA test returned a highly significant relationship (p ~ 0.00) between observed and predicted erosivity values. The R2 statistic (square of the correlation coefficient in Fig. 3a) indicated that the REMM explains ~ 90% of the erosivity variability. MAE (mean absolute error) was equal to 55 MJ mm hm−2 h-1 year−1—which is ~ 3% of the mean erosivity over the period 1994–2015—with the Kling-Gupta Efficiency (KGE) equal to 0.89. The calibrated regression (Fig. 3a, black line, Eq. (1)) shows only negligible departures of data-points from the 1:1 identity line (red line), indicating that the ability of the model (Eq. (1)) to predict actual erosivity in the MedCA is satisfactory. In particular, the regression line has an intercept a = –5.2 (± 95.6) and a slope b = 1.0 (± 0.1) near or equal to the optimum values (a = 0 and b = 1). Though the relatively high standard error of the intercept (± 95.6) indicates the model’s lesser predictive ability for near-zero erosivity values the intercept is not statistically different from zero (Student-t P ~ 1.0). The Nash–Sutcliffe efficiency value obtained in the calibration stage (EF = 0.9) also indicates limited uncertainty in model estimates. The distribution of the residuals approaches the normal distribution shape, indicating skew-free distribution of errors as the data points are mostly aligned along the QQ-plot function (Fig. 3a1). For the Durbin–Watson (DW) statistics (DW = 2.76), there is no indication of serial autocorrelation in the residuals (p = 0.97).Figure 3 Model calibration and validation. (a) Scatterplot of regression model (black line, Eq. (1) and red line of identity) vs. actual rainfall erosivity estimated upon the MedCA over the period 1994–2015, with the inner bounds showing 90% confidence limits (power pink coloured area), and the outer bounds showing 95% prediction limits for new observations (light pink). (a1) related QQ-plot of residuals at calibration stage. (b) Time co-evolution (1981–1993) of actual (black curve) and estimated (red curve) of rainfall erosivity at validation stage (in a,b, r stands for linear correlation coefficient). For the period 1981–1992, the ANOVA p-value less than 0.05 means that a statistically significant relationship between the estimated and actual data is maintained with the independent dataset used for validation. Though the R2 statistic (square of the correlation coefficient in Fig. 3b) indicates that the model only explains ~ 58% of the variability of actual erosivity, the time-variability as a whole was well reproduced (Fig. 3b). The regression parameters a = − 55.4 ± 322.8 and b = 1.1 ± 0.3 indicate that some model overestimation of near-zero erosivity data can occur with the model but the Nash–Sutcliffe efficiency value equal to 0.6 corroborates that also at the validation stage the uncertainties associated with model estimates are not large. The MAE was equal to 84 MJ mm hm−2 h−1 year−1—which is ~ 5% of mean erosivity over the period 1981–1993—while KGE was 0.5. This analysis also indicated that the categorical storm index SSI(GF) was a reliable substitute of storm rainfall. The importance of including this input in the model was confirmed by the values of r, which increased from 0.88 to 0.95 at calibration stage, and from 0.36 to 0.76 at validation stage. Also the mean absolute percentage error (MAPE), which measures the prediction accuracy of a model, indicated an increased accuracy, from 6.3 to 4.3% at calibration stage, and from 9.4 to 7.0%, at validation stage, when including SSI(GF) as input. Historical rainfall erosivity reconstruction Figure 4 shows the areal-mean erosivity evolution over the period 1680–2019 CE, as obtained by means of Eq. (1). The time-series was analysed to find out possible patterns of rainfall aggressiveness and to compare contemporary conditions with historical erosivity patterns (Fig. 4a). The long-term mean value of estimated erosivity data is 1367 ± 319 (SD) MJ mm hm−2 h−1. The application of the Standard Normal Homogeneity Test for the double shift46 (which is useful for dividing long time-series into shorter periods) suggests discontinuities (change-points) in the annual erosivity time-series at a 99% confidence interval in the years 1839 and 1858. The fingerprint of climate change was also detected with the quantile erosivity data with 30-year return period (red curve), with change points in 1838 and 1861, not dissimilar from those obtained with the erosivity data time-series. Other test statistics47–52 detected change points in 1873 and 1903 for the erosivity data, or 1876 and 1906 for the quantiles, which merely support the idea of a long transition period going from the final phase of the Little Ice Age (LIA; ~ 1300–1850 CE53) to the most recent warming. These different statistically-relevant years provide a loose picture of climate-related erosivity variations with changing climate patterns, where the cold conditions of the LIA are still dominating after the end of the Dalton minimum of reduced solar activity (~ 1790–1830 CE) until towards the end of the nineteenth century, but in the process of evolving into an incipient warming that becomes noticeable later in the twentieth century. Following the first detected change point, the quantiles’ values evolve according to a second-order polynomial, whose descending portion passes from 2500 to 1900 MJ mm hm−2 h−1 between the change-point around 1839 and 1950s (white curve on grey band). Our analysis supports a relation between our output (rainfall erosivity) and the dynamics of the North Atlantic Oscillation (NAO)54. The NAO represents a redistribution of air masses between sub-tropical (Azores archipelago, roughly 38°N) and sub-polar latitudes (Iceland, roughly 65°N), and modulates the strength and latitudinal location of the westerly flows, whose major influence on Central Mediterranean precipitation is documented55, 56. Both the positive and negative phases of the NAO are associated with regional changes in large-scale modulations of zonal and meridional heat and moisture transport patterns57. In particular, the positive phases of the NAO reflect lower than normal heights and pressure in the high latitudes of the North Atlantic, and higher than normal heights and pressure in the central North Atlantic and Western Europe (with negative phases reflecting opposite patterns). Strongly positive phases of NAO are usually associated with below-average temperatures and precipitation in southern Europe, while opposite anomalies of temperature and precipitation are generally observed during strongly negative phases. While the NAO shows considerable inter-seasonal and interannual variability, the winter NAO also shows considerable multi-decadal variability58. For instance, the negative phase of the NAO dominated the Atlantic circulation between the mid-1950s to the 1970s. An abrupt transition to recurrent positive phases of the NAO then occurred during the winter 1979–1980 (with the atmosphere still locked into this mode during the winter 1994–1995), followed by a return to a strongly negative phase of the NAO59. Figure 4a shows that phases in the studied period (1680–2019) turning to substantially neutral (with no clear dominance of positive or negative episodes) NAO state (0.12 ± 0.09 standard error, horizontal grey dashed line after the change point) correspond to a decline of the erosivity quantile.Figure 4 Overview of several precipitation patterns over MedCA. (a) Timeline evolution of annual rainfall erosivity reconstructed by Eq. (1) (blue curve) and its quantile with return period (RP) = 30 years (red curve) in a 21-year running window (in MJ mm hm−2 h−1 year−1) and with the related plynomial trend (white curve) after the change-point. Horizontal grey dashed lines are mean values of the North Atlantic Oscillation (NAO) for two periods (reversed axis), calculated from the proxy-based multidecadal winter NAO reconstruction of Trouet et al.100 (orange line). (b) Erosivity density timeline (blue curve) with its quantile with RP = 30 years during the period 1680–2019 (in MJ mm hm−2 h−1 year−1). (c) Precipitation fraction due to very wet days (PF > 95th percentile) for the MedCA (black curve), northern Italy (orange curve) and Sicily (blue curve) over the period 1948–2019 (in mm d−1). (d) Anomalies in the spatial pattern of annual mean rain rates (1991–2019 minus 1961–1990) over the Mediterranean region (from NCEP/NCAR reanalysis data, https://psl.noaa.gov/cgi-bin/data/getpage.pl). We derived that the interannual variability of rainfall erosivity is pronounced either before or after the change-point (Fig. 4). However, the MedCA appears to be subject to more frequent peak values before the break-point (e.g., in 1707 with 2095, in 1757 with 2725, in 1758 with 2760, and in 1812 with 2126 MJ mm hm−2 h−1), a period dominated by the negative mode (− 0.40 ± 0.06 standard error) of the NAO. This is in line with results from previous studies showing that the erosive power of rainfall and extreme precipitation is stronger during cyclonic (low-pressure) conditions (for the western Mediterranean60; and for Montenegro, on the east side of the MedCA61). Hydrological extremes may persist and evolve in unexpected and erratic ways; in fact, the estimated rainfall erosivity indicates a recovery in the low-frequency erosive activity of rains during the past two decades. It is difficult to detect signals of climate change in erosivity extremes associated with torrential rainfall (e.g., hourly-long events), due to the wide spatial variability and unpredictability of these events in long time-series. We addressed this issue by examining the erosivity density (ED), i.e., the erosivity per unit of rain, which is a better indicator of the climate erosive hazard than rainfall erosivity31. The evolution of ED is shown in Fig. 4b, and it is surprisingly characterized by a significant increasing linear trend. The related ED-quantile (RP = 30, red curve) reveals a continuous and strong growth, with a significant long-term linear trend (Mann–Kendall test, p < 0.01). In this case, ED-quantile increases from 1.2, at beginning, to 1.4 MJ hm−2 h−1 year−1, at the end of time-series. Since both erosivity and erosive density appear subject to a complex evolution over the most recent (warmest) decades. We inspected the period 1948–2019 in more detail by arranging the precipitation fraction (PF) due to very wet days (PF > 95th percentile) from the NCEP/NCAR Reanalysis data. We found that the PF > 95th percentile improves the accuracy of heavy precipitation estimates of very wet days to total precipitation from the probability distribution of daily precipitation than from the raw data62. In this way, we detected a temporal development of the P > 95 with an intensifying trend during the recent warming phase, 1986–2019, at both regional and sub-regional scales, confirming how the landscape, during recent decades, have been recurrently subjected to gradually increasing hydrological stress (Fig. 4c). Black curve in Fig. 4c illustrates the evolution of FP > 95 over the MedCA, while blue and orange curves represent the FP > 95 for northern and insular Italy (Sicily), respectively. The anomaly map furthermore suggests an enhanced hazard associated with complex and more intense rainfall events over the MedCA during recent decades (Fig. 4d). It is characterised by a strong convective component, especially in the autumn season63. As a consequence, disaster-affected areas in the MedCA have become more exposed to climate hazard conditions because rainfall aggressiveness has apparently become increasingly changeable and unpredictable at small scales64. According to Cislaghi et al.65 and Pavan et al.66, the frequency of occurrence of daily precipitation has decreased over Italy, but short-duration episodes (i.e., from 1 to 3 h) have instead enhanced the torrential character of seasonal rains. Colarieti Tosti67 also reported that in the coming decades the polar vortex would likely go through a phase of expansion towards the southern latitudes, with consequent exacerbation of the hydrological cycle in the Mediterranean. An increasing frequency of extreme precipitation events is also expected over hazard-exposed landscapes over much of Europe towards the end of the twenty-first century9, 68. Climate model simulations with the new generation of Coordinated Downscaling Experiment over Europe (EURO-CORDEX69) reveal an increasing trend towards a higher frequency of hydrologic 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 for most European countries—including southern Europe—no significant change in annual precipitation, but at the same time an increase in maximum daily precipitation70. These results seem to reflect the paradoxical increase of Mediterranean extreme rainfall in spite of decrease in total values, as claimed by Alpert et al.71, Paxian et al.72 and Caloiero et al.73 but questioned by Mariani and Parisi74 for the whole basin and the western sub-basin. Other researchers75 also indicated no significant trends in the observed data, testifying a substantial stability of the temporal and spatial behaviour of heavy rain events over the course of the twentieth century. These contrasting results highlight that different analyses can change perception of the type of hazard associated with hydrological processes when using dissimilar metrics. With this study, we advocate the use of metrics that not only reflect the climate forcing component reproduced in the prevailing storm aggressiveness (erosivity), but in addition the damaging hydrologic hazard. The erosivity density offers this opportunity because regions with high ED values are exposed to a risk of flooding (and even water scarcity) because of their infrequent, but very intense and erosive rainstorms76. While for the twentieth century, individual time-series of erosivity data are available at some Italian sites, a reconstruction was lacking that would enable a comparison between recent and historical erosivity covering particular climatic periods like the LIA. For that, we have developed and assessed a historical model (traced back to the seventeenth century), which can explain the erosivity data made available from previous studies. In this way, we have framed the historical trend of erosivity (with a view to exploring the transition from the LIA to the modern warming) rather than reconstructing with an arguable level of confidence the erosivity in each single year. Our work reveals that improved knowledge of the temporal changes in erosive precipitation is critical since these changes indicate different degrees of landscape damages and are of great importance for the implementation of environment conservation and management plans under a changing climate77. The occurrence of erosivity events is expected to be significantly altered under a warming climate. Since southern Europe landscapes tend to react rapidly to such changes, understanding extreme precipitation and the associated erosive conditions is vital to better comprehend and anticipate environmental changes in this area. Meeting the demand for reliable projections of extreme, short-duration rainfall is challenging because of our poor understanding of past hydroclimatic dynamics and mechanisms, which limits our ability to generate plausible predictions about future states11. As we have shown in this study, precipitation and storm indicators provide new opportunities to develop long records of past erosivity data. We improved the understanding of hydrological extremes, and the inter-annual variability of erosive rainfall and erosivity density, employing long continuous time-series covering 1680–2019 CE, developed for the MedCA from available instrumental and documentary climate records. The integrated approach used to reconstruct and assess yearly erosivity data allowed detecting signals of present-day climate change. For the entire assessed period 1680–2019 CE, the main conclusions can be summarized as follows:Annual rainfall erosivity data across the MedCA from 1680 onwards show alternate stormy and quieter periods, with a change point towards the end of the Little Ice Age (1839) and different trends before (1680–1838 CE) and after (1839–2019 CE) that break point. Conversely, erosivity density has undergone an increasing trend since the beginning of the series. The generally cold period 1680–1838 CE appears characterized by some peaks of stormy years, followed by a quieter phase, marked by a transition to decreasing storminess over the dominant mode of neutral NAO. During recent decades, however, the precipitation fraction due to very wet days (> 95th percentile) and rainfall erosivity indicate a recovery in the erosive activity from more intensive rainfall, concurrent with an increased variability of erosive density. The way how in the recent decades Mediterranean cyclones have been producing trends of erosivity from rainfall extremes remains elusive24, 25. Whether cyclone occurrences decreased over the Mediterranean, a rising erosivity associated with short-term and very extreme rain events could also be due to a strong increase of extreme related convective-precipitation during recent warming. Berg et al.78 showed that convective precipitation is more responsive to temperature increases than frontal passages, and increasingly dominates rain events of extreme intsensity. Chains of damaging events would likely become more common with increasing climate warming79, though with distinct patterns of change in small domains without the emergence of spatially and temporally homogeneous trends80. The present analysis, performed at an annual temporal scale, may mask important variations manifesting at even finer time-scales. The methodology applied mostly addresses inter-annual and inter-decadal time-scales, and does not capture daily to seasonal changes, which may impact on hydrology and lead to damages due to losses of land at different scales. For the MedCA, large erosive rain events were observed to occur especially in summer in continental areas, and in autumn along the coasts and near-to-coast reliefs81. It was also shown82 that more frequent extreme events in autumn do not cause seasonal rainfall totals to deviate from the historical range of climate variation. They rather tend to generate more disproportion between dry and wet periods, which could bring soil loss to higher rates83. The regional analysis may also obscure sub-regional trends84. However, only a few sub-regional studies extend back erosivity data over such a long period as in our study. Erosivity time-series reconstructed in two Mediterranean fluvial basins, the Calore river basin (southern Italy)85 and in the Po river basin (northern Italy)86, showed similar increasing trends since the end of the LIA. Diodato et al.87 also showed an upward trend in the Swiss plateau (north of the Mediterranean central area). If this rainfall regime would continue, it could result in an ever-increasing erosive hazard affecting Mediterranean lands in a more erratic fashion. This implies that a future increase in extreme precipitation events might have severe consequences regarding soil erosion, flash-flood risk, and various types of ecological disruptances in the Mediterranean region88. Furthermore, this underlines the need for high-resolution climate model projections to reliably predict erosive precipitation. In conclusion, by providing the first quantitative assessment of the long-term dynamics of storm-related extremes in the MedCA, this work highlights the promise of this approach and provides impetus for further refinement of erosive rainfall reconstructions. The time-series of annual erosivity and erosive density presented here are arguably of great importance for climatic analyses dealing with climate variability at multidecadal to centennial time-scales. It is on precisely these time-scales that the anthropogenic-forced changes are most likely to be superimposed on natural (forced and unforced) climate variability. The erosivity data presented here can be used to study the variability and the extremes of climate occurring at a sub-regional scale, which can be compared to model-based simulations of natural and forced (external and internal) variability for the past centuries. This study also indicates the importance of erosivity density for representing and predicting future hydrological changes across the Mediterranean region. The limited availability of basin-wide consistent erosivity data constrained the ability of our study to comprehensively model erosivity for the entire region. Furthermore, with this limitation in mind, these results provide essential information to shed light on the processes governing long-term storm-related phenomena. This implies that environmental management can use data from long historical time-series as a basis for increasing societal resilience, informing decision making and directing new monitoring/modelling efforts. In particular, they can raise the awareness of policymakers by pointing to the urgent need of minimizing and controlling the land degradation in the Mediterranean basin. Finally, our study emphasises the importance of separating low-frequency trends from high-frequency extreme events that, as revealed here, may show entirely opposite trends. The inability of most early instrumental time-series, and climate reconstructions, to capture high-frequency extreme events have clearly limited the possibility to fingerprint anthropogenic climate change against the full range of long-term natural variability across time-scales. Whereas the mean state, and rate of change, for a particular climate variability may still be within the range of natural variability, the amplitude or return rate of certain climate extremes may not. Methods Study area The MedCA has a roughly rectangular shape centred on Italy, with a preferential north–south orientation (36°–47° N, 5°–18° E). It extends over large portions of southern Europe, in the transition zone between Central Europe and the southern Mediterranean rim. The Alpine sector forms an imposing mountain range to the north, with a length of ~570 km and a width varying from 25 to 110 km. The Apennine mountain chain runs over ~ 1000 km along the Italian peninsula. The Italy islands of Sardinia and Sicily close this geographical variety. The climate of the MedCA is highly diversified, given its geographical position and the varied morphology of the sectors composing it. The meteorological conditions largely depend on the fronts arising as a result of cold polar air meeting warm tropical air. In particular, different situations are produced depending 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. The MedCA is subject to pronounced rainfall erosivity in the context of southern Europe. (R)USLE-based actual rainfall erosivity data Annual (R)USLE-based R-factor data (1981–2015) were derived from Diodato89 and Acquaotta et al.90, with updating from the network of SCIA (http://www.scia.isprambiente.it/wwwrootscia/scia.html)—National System for the Collection and Elaboration of Climatological Data91, for a total of 37 stations (Fig. 6). For the period 1994–2015, (R)USLE-based data were available for the full set of meteorological stations (Fig. 6, both black and red numbers), because the digital National Agrometeorological Network (https://tinyurl.com/h4juzuv) came into operation in 1994. Prior to 1994, only some data were available for the Piedmont Region of Italy (north-western part of the MedCA), with only few sparse stations providing data since 1981. The most accurate dataset, 1994–2015, was used for model calibration, while the least accurate dataset, 1981–1993, was used for model validation (Fig. 5, red numbers).Figure 5 Spatial distribution of digital rain gauges used for calibration (black numbers) and validation (red numbers) across the Mediterranean Central Area. Numerical and categorical inputs The reconstruction of annually-resolved rainfall erosivity for the MedCA was based on seasonal precipitation data covering most of Europe (30° W–40° E/30° N–71° N) from 1500–2000 on 0.5° by 0.5° grid resolution92. The GPCC (Global Precipitation Climatology Centre) provides an extension of the seasonal precipitation dataset until 2019 (retrieved from https://tinyurl.com/wlnvznq). GPCC products are known to outperform other similar products (e.g. CRU, ERA‐Interim, ERA‐40) in Mediterranean areas93, 94. A categorical variable, the annual severity storm index sum (SSI) from Diodato et al.26 was used as a proxy to overcome the lack of reliable information about historical rain intensity. SSI was derived from several written 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). The SSI time-series was smoothed by applying the low-pass Gaussian filtering technique. Rainfall erosivity mediterranean model For the historical reconstruction of annual rainfall erosivity (MJ mm hm−2 h−1 year−1), we developed a regression model, hereafter referred to as the Rainfall Erosivity Mediterranean Model (REMM), which uses seasonal precipitation data (mm) and the Gaussian-filtered Annual Storm Severity Index (SSI(GF) as inputs. The non-linear dependence of rainfall erosivity on rain intensity supports the adoption of nonlinear indicators within a parsimonious approach comparable to the (R)USLE approach4. The non-linear model takes the following form: 1 REMM=A·SDMax(Ps)·PSum·PAut+B·α+ASSIS(GF)k·PWin+PSpr+C where A (MJ mm−1 hm−2 h−1 year−1) and B (MJ hm−2 h−1 year−1) are scale parameters converting the result of multiplications into the output unit, and C (MJ mm hm−2 h−1 year−1) is a shift parameter estimating rainfall erosivity when the seasonal precipitation (P, mm) inputs (Sum: summer, Aut: autumn, Win: winter, Spr: spring) are equal to zero. The intensity of the product between summer and autumn precipitation (PSum·PAut) is modulated by a modified version of the variation coefficient, where SD is the inter-seasonal standard deviation and Max(Ps) is the maximum seasonal precipitation per year (mm). The product (PSum·PAut) was considered here because the power of rainfall exerts the maximum of erosive forces in summer and autumn (Fig. 6).Figure 6 Seasonal regime of storm cells with the related scheme of variability of the monthly storm erosivity strength over the MedCA. Severe convective erosive storms occur between June and November in Europe32. For this region, the occurrence of severe thunderstorms is indeed intimately associated with the convective environmental conditions. The frequencies of flood and flash-flood events shows that the precipitation in spring likely leads to a different regime of hydrological extremes compared to June–October, when flash-floods are more frequent32. This means that the period June–November delineates a key time window in estimating rainfall–runoff erosivity in Europe. The sum of winter and spring rainfalls (PWin+PSpr) is instead modulated by the storm index SSI(GF) within the φ term α+ASSIS(GF)k. The central Mediterranean, although located to the south of the main Atlantic storm track that more directly affects northwestern Europe, is quite frequently subject to sudden events of extreme and adverse weather, which cannot be accounted for by only seasonal precipitation totals26. The additive component PWin+PSpr supports the assumption of a softer dependence of erosivity on precipitation intensity. The concept of the model is summarised in Fig. 7 following Waldam95 and Diodato and Bellocchi96. A separation of advective and convective events is considered important, assuming that convective precipitation (brief and intense) has higher dynamics and variability than the usually more static advective events. Seasonal differences in precipitation patterns are governed either by convective processes, i.e. higher dynamics more common in summer (when rain-splash dominate), or a variable mix of convective and advective rains more frequent in autumn (when diffuse overland flow dominates). Then, winter and spring rains (of long duration and low intensity), usually originating from broad mid-latitude frontal activity, carry high volumes of rainwater thanks to orographically enhanced stratiform precipitation causing large-scale hydrological processes like flooding97. High percentiles are used to differentiate more clearly such rainfall characteristics These processes are captured in the Rainfall Erosivity Mediterranean Model—Eq. (1)—by the storm-severity index (ASSIS(GF)). Highly intensive convective (or of a mixed advective-convective nature) rain events result instead in the splash-erosivity process, as interpreted by the variation coefficient SDMax(Ps). In this case, raindrops with high kinetic energy may cause local downpours.Figure 7 Monthly precipitation (98th percentiles) and frequency of floods and flash-floods (bars) driven by dominating frontal, convective or advective-convective rainfall events over the Mediterranean basin. The number of monthly floods (blue bars) and flash-floods (black bars) is from Guzzetti et al.101, Gaume et al.102 and Diodato et al.103 for the period 1950–2006. For the same period, monthly percentiles (prc) were extracted from the CRU Global Climate Dataset104. The terms of Eq. (1) are reported. Trend assessment of extreme values The quantile approach was applied to identify step trends at any EDX scheme with assigned return periods. The return period (T) for the seasonal rainfall erosivity density falling above the jth quantile was ranked using the lognormal distribution: 5 QEDXt-wT=expμEDXt-w∗+uT·σEDXt-w∗ where QEDXt-wT is the jth rainfall erosivity density (EDX) quantile of the log-normal distribution with assigned return period; μEDXt-w∗ and σEDXt-w∗ are the mean and the standard deviation of the variable EDX t* = ln(EDX t). The subscript t-w indicates the computation of the generic variable at the t-time over a 22-year moving window (w); uT is the log-normal dimensionless coefficient equal to 1.27 for T = 30 years. The moving window of ~ 22-years reflects the periodicity of climate phenomena, which proved effective to show long-term trends by smoothing out the secular variation and the more volatile year-to-year changes70. Model calibration and assessment To assess the model, statistical analyses were performed with STATGRAPHICS (http://www.duke.edu/~rnau/sgwin5.pdf), with the graphical support of WESSA (https://www.wessa.net) and CurveExpert routines (https://www.curveexpert.net). The parameters of Eq. 1 were calibrated against actual rainfall erosivity data according to statistical criteria. The first condition was to minimize the distance between modelled and actual erosivity data, by minimizing the Mean Absolute Error (optimum, 0 ≤ MAE < ∞, MJ mm hm-2 h-1 yr-1). Complementary to the MAE, the MAPE (mean absolute percent error) offers the advantage of being scale-independent and intuitive (e.g. the prediction model is considered reasonable with a MAPE below 30% and very accurate with a MAPE less than 10%). The second condition is to maximise the determination coefficient (0 ≤ R2 ≤ 1, optimum) that is the variance explained by the model. The third conditions approximates the unit slope of the straight line that would minimise the bias of the linear regression actual versus modelled data (b = 1, optimum). In addition, the Kling-Gupta index (–∞ < KGE ≤ 1) was used as efficiency measure, with KGE > –0.41 indicating that a model improves upon the means of observations as a benchmark predictor. The Nash–Sutcliffe efficiency (-∞ < EF ≤ 1, optimum) was also calculated as an uncertainty indicator of the model performance because greater values than 0.6 indicate limited model uncertainty, likely associated with narrow parameter uncertainty98. To select the set of important covariates for the parsimonious model for estimating actual erosivity data, we iteratively added in predictors, one-at-a-time until modelling solutions with small MAE and large R2 values were obtained. Then, for the final selection, the third criterion—|b-1|= min—was additionally involved. Each predictor was repositioned over > 50 iterations until convergence was achieved99. The Durbin-Watson statistic was performed to test for auto-correlated residuals because large temporal dependence may induce spurious correlations. ANOVA p-values were used to present the statistical significance of the regression between estimates and the actual data. Supplementary Information Supplementary Information Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-020-78857-3. Acknowledgements N. D. and G. B. performed this research as an investigator-driven study without financial support. F. C. L. was supported by the Swedish Research Council (Vetenskapsrådet, grant no. 2018-01272), and conducted the work with this article as a Pro Futura Scientia XIII Fellow funded by the Swedish Collegium for Advanced Study through Riksbankens Jubileumsfond. The publication cost was covered by Stockholm University. Author contributions N.D. and G.B. developed the original research design and collected and analysed the historical documentary data. N.D., F.C.L. and G.B. wrote the article together and made the interpretations together. All authors reviewed the final manuscript. Funding Open Access funding provided by Stockholm University. Data availability All data used in this study are freely available. Spatial patterns of mean annual rainfall erosivity over the European region (Fig. 1c) are freely available from the ESDAC (European Soil Data Centre) dataset at https://esdac.jrc.ec.europa.eu/content/global-rainfall-erosivity. Anomalies in the spatial pattern of annual mean rain rates over the Mediterranean region (Fig. 4d) are obtained from NCEP (US National Center for Environmental Prediction)/NCAR (US National Center for Atmospheric Research) reanalysis data at https://psl.noaa.gov/cgi-bin/data/getpage.pl. The proxy-based multi-decadal winter NAO reconstruction (Fig. 4a) is by Climate Explorer Climate Change Atlas of the Dutch Royal Netherlands Meteorological Institute (KNMI) at http://climexp.knmi.nl. Updated (R)USLE-based erosivity data (Fig. 5) were derived from the system for the collection and elaboration of climatological data (SCIA) of the Italian National Institute for Environmental Protection and Research (ISPRA) at http://www.scia.isprambiente.it/wwwrootscia/scia.html, completed with the data provided from the digital National Agrometeorological Network accessible through https://tinyurl.com/h4juzuv. Updated seasonal precipitation datasets were retrieved from the Global Precipitation Climatology Centre (GPCC) through https://tinyurl.com/wlnvznq. Also, the full set of raw data and the equations that support the findings of this study (erosivity model, time-series reconstruction and precipitation fraction data, seasonal precipitation and storm-severity index inputs), are available in the supplementary Table S1. Competing interests The authors declare no competing interests. ==== Refs References 1. Wischmeier, W. H. & Smith, D. D. Predicting rainfall erosion losses: A guide to conservation planning (Washington, DC: U.S. Department of Agriculture, Agriculture Handbook No. 537, 1978). 2. Brown LC Foster GR Storm erosivity using idealized intensity distributions Trans. ASABE 1987 30 379 386 10.13031/2013.31957 3. Reinard KG Freimund JR Using monthly precipitation data to estimate the R-factor in the revised USLE J. Hydrol. 1994 157 287 306 10.1016/0022-1694(94)90110-4 4. 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) (Washington, DC: USDA-ARS Agriculture Handbook No. 703, 1997). 5. Mondal A Khare D Kundu S Change in rainfall erosivity in the past and future due to climate change in the central part of India Int. Soil Water Conserv. Res. 2016 4 186 194 10.1016/j.iswcr.2016.08.004 6. Toreti A Cronie O Zampieri M Concurrent climate extremes in the key wheat producing regions of the world Sci. Rep. 2019 9 5493 10.1038/s41598-019-41932-5 30940858 7. Jongman B Effective adaptation to rising flood risk Nat. Commun. 2018 9 1986 10.1038/s41467-018-04396-1 29844334 8. Hoeppe P Trends in weather related disasters—Consequences for insurers and society Weather Clim. Extrem. 2016 11 70 79 10.1016/j.wace.2015.10.002 9. Myhre G Frequency of extreme precipitation increases extensively with event rareness under global warming Sci. Rep. 2019 9 16063 10.1038/s41598-019-52277-4 31690736 10. Boudet H Giordono L Zanocco C Satein H Whitley H Event attribution and partisanship shape local discussion of climate change after extreme weather Nat. Clim. Change 2020 10 69 76 10.1038/s41558-019-0641-3 11. Zhang X Zwiers FW Li G Wan H Cannon AJ Complexity in estimating past and future extreme short-duration rainfall Nat. Geosci. 2017 10 255 259 10.1038/ngeo2911 12. Blenkinsop S The INTENSE project: using observations and models to understand the past, present and future of sub-daily rainfall extremes Adv. Sci. Res. 2018 15 117 126 10.5194/asr-15-117-2018 13. Sillmann J Sippel S Russo S Climate Extremes and Their Implications for Impact and Risk Assessment 2019 Amsterdam Elsevier 14. Ljungqvist FC Krusic PJ Sundqvist HS Zorita E Brattström G Frank D Northern Hemisphere hydroclimate variability over the past twelve centuries Nature 2016 532 94 98 10.1038/nature17418 27078569 15. Xie S-P Towards predictive understanding of regional climate change Nat. Clim. Change 2015 5 921 930 10.1038/nclimate2689 16. Schewe J State-of-the-art global models underestimate impacts from climate extremes Nat. Commun. 2019 10 1005 10.1038/s41467-019-08745-6 30824763 17. D’Asaro F D’Agostino L Bagarello V Assessing changes in rainfall erosivity in sicily during the twentieth century Hydrol. Process. 2007 21 2862 2871 10.1002/hyp.6502 18. Toy TJ Foster GR Renard KG Soil Erosion; Prediction, Measurement, and Control 2002 New York John Wiley & Sons Inc 19. Reimann L Vafeidis AT Sally Brown A Hinkel J Tol RSJ Mediterranean UNESCO World Heritage at risk from coastal flooding and erosion due to sea-level rise Nat. Commun. 2018 9 4161 10.1038/s41467-018-06645-9 30327459 20. Schmidt S Alewell C Panagos P Meusburger K Regionalization of monthly rainfall erosivity patterns in Switzerland Hydrol. Earth Syst. Sc. 2016 20 4359 4373 10.5194/hess-20-4359-2016 21. Zittis G Hadjinicolaou P Klangidou M Proestos Y Lelieveld J A multi-model, multi-scenario, and multi-domain analysis of regional climate projections for the Mediterranean Reg. Environ. Change 2019 19 2621 2635 10.1007/s10113-019-01565-w 22. Lenderink G Fowler HJ Understanding rainfall extremes Nat. Clim. Change 2017 7 391 393 10.1038/nclimate3305 23. Diodato N Estimating RUSLE’s rainfall factor in the part of Italy with a Mediterranean rainfall regime Hydrol. Earth Syst. Sci. 2004 8 103 107 10.5194/hess-8-103-2004 24. Lionello P Objective climatology of cyclones in the Mediterranean region: A consensus view among methods with different system identification and tracking criteria Tellus A Dyn. Meteorol. Oceanogr. 2016 68 29391 10.3402/tellusa.v68.29391 25. Camuffo D Jones PD Improved understanding of past climatic variability from early daily European instrumental sources Clim. Change 2002 53 1 392 10.1023/A:1014902904197 26. Diodato N Ljungqvist FC Bellocchi G A millennium-long reconstruction of damaging hydrological events across Italy Sci. Rep. 2019 9 9963 10.1038/s41598-019-46207-7 31292466 27. Pichard G Arnaud-Fassetta G Moron V Roucaute E Hydroclimatology of the Lower Rhône Valley: historical flood reconstruction (AD 1300–2000) based on documentary and instrumental sources Hydrol. Sci. J. 2017 62 1772 1795 10.1080/02626667.2017.1349314 28. Benito G Recurring flood distribution patterns related to short-term Holocene climatic variability Sci. Rep. 2015 5 16398 10.1038/srep16398 26549043 29. Allan R Toward integrated historical climate research: The example of atmospheric circulation reconstructions over the earth Wiley Interdiscip. Rev. Clim. Change 2016 7 164 174 10.1002/wcc.379 30. Diodato N Borrelli P Panagos P Bellocchi G Bertolin C Communicating hydrological hazard-prone areas in Italy with geospatial probability maps Front. Environ. Sci. 2019 7 193 10.3389/fenvs.2019.00193 31. Glur L Frequent floods in the European Alps coincide with cooler periods of the past 2500 years Sci. Rep. 2013 3 2770 10.1038/srep02770 24067733 32. Diodato N Bellocchi G Decadal modelling of rainfall–runoff erosivity in the Euro-Mediterranean region using extreme precipitation indices Glob. Planet. Change 2012 86–87 79 91 10.1016/j.gloplacha.2012.02.002 33. Rogers JC North Atlantic storm track variability and its association to the North Atlantic oscillation and climate variability of Northern Europe J. Clim. 1997 10 1635 1647 10.1175/1520-0442(1997)010<1635:NASTVA>2.0.CO;2 34. Xoplaki E Modelling climate and societal resilience in the Eastern Mediterranean in the last millennium Hum. Ecol. 2018 46 363 379 10.1007/s10745-018-9995-9 35. Dayan U Nissen K Ulbrich U Atmospheric conditions inducing extreme precipitation over the eastern and western Mediterranean Nat. Hazard. Earth Syst. 2015 15 2525 2544 10.5194/nhess-15-2525-2015 36. Hofstätter M Chimani B Lexer A Blösch A new classification scheme of European cyclone tracks with relevance to precipitation Water Resour. Res. 2016 52 7086 7104 10.1002/2016WR019146 37. Diodato N Ceccarelli M Bellocchi G Decadal and century-long changes in the reconstruction of erosive rainfall anomalies at a Mediterranean fluvial basin Earth Surf. Proc. Land. 2008 33 2078 2093 10.1002/esp.1656 38. Poirier C Poitevin C Chaumillon E Comparison of estuarine sediment record with modelled rates of sediment supply from a western European catchment since 1500 C. R. Geosci. 2016 348 479 488 10.1016/j.crte.2015.02.009 39. Pinto, J. C., Klawa, M., Ulbrich, U., Rudari, R. & Speth, P. Extreme precipitation events over northwest Italy and their relationship with tropical-extratropical interactions over the Atlantic. In: Mediterranean storms (eds Deidda, R., Mugnai, A. & Siccardi, F.) 1–6 (2001). 40. Wainwright, J. Weathering, soils and slope processes. In: The physical geography of the Mediterranean (ed Woodward, J. C.) 169–202 (2009). 41. López-Vicente M Navas A Machín J Identifying erosive periods by using RUSLE factors in mountain fields of the Central Spanish Pyrenees Hydrol. Earth Syst. Sc. 2008 12 523 535 10.5194/hess-12-523-2008 42. Borrelli P Diodato N Panagos P Rainfall erosivity in Italy: A national scale spatio-temporal assessment Int. J. Digit. Earth 2016 9 835 850 10.1080/17538947.2016.1148203 43. Hernando D Romana MG Estimate of the (R)USLE rainfall erosivity factor from monthly precipitation data in mainland Spain J. Iber. Geol. 2016 42 113 124 10.5209/rev_JIGE.2016.v42.n1.49120 44. Ballabio C Mapping monthly rainfall erosivity in Europe Sci. Total Environ. 2017 579 1298 1315 10.1016/j.scitotenv.2016.11.123 27913025 45. Diodato N Borrelli P Fiener P Bellocchi G Romano N Discovering historical rainfall erosivity with a parsimonious approach: A case study in Western Germany J. Hydrol. 2017 544 1 9 10.1016/j.jhydrol.2016.11.023 46. Wanner H North Atlantic oscillation—concepts and studies Surv. Geophys. 2001 22 321 382 10.1023/A:1014217317898 47. Alexandersson H Moberg A Homogenization of Swedish temperature data. Part I: homogeneity test for linear trends Int. J. Climatol. 1997 17 25 34 10.1002/(SICI)1097-0088(199701)17:1<25::AID-JOC103>3.0.CO;2-J 48. Pettitt AN A non-parametric approach to the change-point problem Appl. Stat. 1979 28 126 135 10.2307/2346729 49. Buishand TA Some methods for testing the homogeneity of rainfall records J. Hydrol. 1982 58 11 27 10.1016/0022-1694(82)90066-X 50. Alexandersson H A homogeneity test applied to precipitation data J. Climatol. 1986 6 661 675 10.1002/joc.3370060607 51. Worsley KJ Confidence regions for a change-point in a sequence of exponential family random variables Biometrika 1986 71 91 104 10.1093/biomet/73.1.91 52. Wang XL Wen QH Wu Y Penalized maximal t test for detecting undocumented mean change in climate data series J. Appl. Meteorolog. Clim. 2007 46 916 931 10.1175/JAM2504.1 53. Miller GH Abrupt onset of the Little Ice Age triggered by volcanism and sustained by sea-ice/ocean feedbacks Geophys. Res. Lett. 2012 39 L02708 10.1029/2011GL050168 54. Wagner S Zorita E The influence of volcanic, solar and CO2 forcing on the temperatures in the Dalton minimum 1790–1830: A model study Clim. Dyn. 2005 25 205 218 10.1007/s00382-005-0029-0 55. Trigo RM Osborn TJ Corte-Real JM The North Atlantic oscillation influence on Europe: Climate impacts and associated physical mechanisms Clim. Res. 2002 20 9 17 10.3354/cr020009 56. Pasini A Langone R Attribution of precipitation changes on a regional scale by neural network modeling: A case study Water 2010 2 321 332 10.3390/w2030321 57. Hurrell JW Decadal trends in the North Atlantic oscillation: Regional temperatures and precipitation Science 1995 269 676 679 10.1126/science.269.5224.676 17758812 58. Chelliah M Bell GD Tropical multidecadal and interannual climate variability in the NCEP–NCAR reanalysis J. Clim. 2004 17 1777 1803 10.1175/1520-0442(2004)017<1777:TMAICV>2.0.CO;2 59. Halpert MS Bell GD Climate assessment for 1996 Bull. Am. Meteor. Soc. 1997 78 S1 S49 10.1175/1520-0477-78.5s.S1 60. Angulo-Martínez M Beguería S Do atmospheric teleconnection patterns influence rainfall erosivity? A study of NAO, MO and WeMO in NE Spain, 1955–2006 J. Hydrol. 2012 450–451 168 179 10.1016/j.jhydrol.2012.04.063 61. Milošević DD Variability of seasonal and annual precipitation in Slovenia and its correlation with large-scale atmospheric circulation Open Geosci. 2016 8 593 605 10.1515/geo-2016-0041 62. Zolina O Simmer C Belyaev K Kapala A Gulev S Improving estimates of heavy and extreme precipitation using daily records from European rain gauges J. Hydrometeorol. 2009 10 701 716 10.1175/2008JHM1055.1 63. Diodato N Bellocchi G Chirico GB Romano N How the aggressiveness of rainfalls in the Mediterranean lands is enhanced by climate change Clim. Change 2011 108 591 599 10.1007/s10584-011-0216-4 64. Diodato, N., Gómara, I. & Bellocchi, G. Recalling the past of erosive rainfall hazard: A tandem with the future projections? Clim. Change (submitted). 65. Cislaghi M De Michele C Ghezzi A Rosso R Statistical assessment of trends and oscillations in rainfall dynamics: Analysis of long daily Italian series Atmos. Res. 2005 77 188 202 10.1016/j.atmosres.2004.12.014 66. Pavan V Tomozeiu R Cacciamani C Di Lorenzo M Daily precipitation observations over Emilia-Romagna: Mean values and extremes Int. J. Climatol. 2008 28 2065 2079 10.1002/joc.1694 67. Colarieti Tosti, C. Il clima del futuro? La chiave è nel passato (https://tinyurl.com/u5hp3w6, 2014). (in Italian) 68. Rädler AT Groenemeijer PH Faust E Sausen R Púčik T Frequency of severe thunderstorms across Europe expected to increase in the 21st century due to rising instability NPJ Clim. Atmos. Sci. 2019 2 30 10.1038/s41612-019-0083-7 69. Jacob D EURO-CORDEX: New high-resolution climate change projections for European impact research Reg. Environ. Change 2014 14 563 578 10.1007/s10113-013-0499-2 70. Alfieri L Burek P Feyen L Forzieri G Global warming increases the frequency of river floods in Europe Hydrol. Earth Syst. Sci. 2015 19 2247 2260 10.5194/hess-19-2247-2015 71. Alpert P The paradoxical increase of Mediterranean extreme daily rainfall in spite of decrease in total values Geophys. Res. Lett. 2002 29 135 154 10.1029/2001GL013554 72. Paxian A Hertig E Seubert S Vogt G Jacobeit J Paeth H Present-day and future mediterranean precipitation extremes assessed by different statistical approaches Clim. Dyn. 2015 44 845 860 10.1007/s00382-014-2428-6 73. Caloiero T The trend of monthly temperature and daily extreme temperature during 1951–2012 in New Zealand Theor. Appl. Climatol. 2017 129 111 127 10.1007/s00704-016-1764-3 74. Mariani, L. & Parisi, S. G. Extreme rainfalls in the Mediterranean area. In: Storminess and environmental change—Climate forcing and responses in the Mediterranean region (eds Diodato, N. & Bellocchi, G.) 17–37 (2014). 75. De Luca DL Galasso L Stationary and non-stationary frameworks for extreme rainfall time series in southern Italy Water 2018 10 1477 10.3390/w10101477 76. Panagos P Ballabio C Borrelli P Meusburger K Spatio-temporal analysis of rainfall erosivity and erosivity density in Greece CATENA 2016 137 161 172 10.1016/j.catena.2015.09.015 77. Liu S Spatial-temporal changes of rainfall erosivity in the loess plateau, China: Changing patterns, causes and implications CATENA 2018 166 279 289 10.1016/j.catena.2018.04.015 78. Berg P Moseley C Haerter JO Strong increase in convective precipitation in response to higher temperatures Nat. Geosci. 2013 6 181 185 10.1038/ngeo1731 79. AghaKouchak A How do natural hazards cascade to cause disasters? Nature 2018 561 458 460 10.1038/d41586-018-06783-6 30250152 80. Libertino A Ganora D Claps P Evidence for increasing rainfall extremes remains elusive at large spatial scales: The case of Italy Geophys. Res. Lett. 2019 46 7437 7446 10.1029/2019GL083371 81. Diodato N Bellocchi G MedREM, a rainfall erosivity model for the Mediterranean region J. Hydrol. 2010 387 119 127 10.1016/j.jhydrol.2010.04.003 82. Diodato N Bellocchi G Assessing and modelling changes in rainfall erosivity at different climate scales Earth Surf. Process. Land. 2009 34 969 980 10.1002/esp.1784 83. Diodato N Bellocchi G Drought stress patterns in Italy using agro-climatic indicators Clim. Res. 2008 36 53 63 10.3354/cr00726 84. Zittis G Observed rainfall trends and precipitation uncertainty in the vicinity of the Mediterranean, Middle East and North Africa Theor. Appl. Climatol. 2018 134 1207 1230 10.1007/s00704-017-2333-0 85. Diodato N Ceccarelli M Bellocchi G Decadal and century-long changes in the reconstruction of erosive rainfall anomalies at a Mediterranean fluvial basin Earth Surf. Process. Landf. 2018 33 2078 2093 10.1002/esp.1656 86. Diodato, N., Ljungqvist, F. C. & Bellocchi, G. Historical predictability of rainfall erosivity – a new reconstruction for monitoring extremes over Northern Italy (1500–2019 CE). NPJ Clim. Atmos. Sci. (in press). 87. Diodato, N., Bellocchi, G., Meusburger, K. & Buttafuoco, G. Modelling long-term storm erosivity time-series: a case study in the Western Swiss Plateau. In: Storminess and environmental change—Climate forcing and responses in the Mediterranean region (eds Diodato, N. & Bellocchi, G.) 149–164 (2014). 88. Borrelli P An assessment of the global impact of 21st century land use change on soil erosion Nat. Commun. 2017 8 2013 10.1038/s41467-017-02142-7 29222506 89. Diodato N Predicting RUSLE (Revised Universal Soil Loss Equation) monthly erosivity index from readily available rainfall data in Mediterranean area Environmentalist 2005 25 63 70 90. Acquaotta F Baronetti A Bentivenga M Fratianni S Piccarreta M Estimation of rainfall erosivity in Piedmont (Northwestern Italy) by using 10-minute fixed-interval rainfall data Idojaras 2019 123 1 18 91. Desiato F Fioravanti G Fraschetti P Perconti W Toreti A Climate indicators for Italy: Calculation and dissemination Adv. Sci. Res. 2011 6 147 150 10.5194/asr-6-147-2011 92. Pauling A Luterbacher J Casty C Wanner H Five hundred years of gridded high-resolution precipitation reconstructions over Europe and the connection to large-scale circulation Clim. Dyn. 2006 26 387 405 10.1007/s00382-005-0090-8 93. Belo-Pereira M Dutra E Viterbo P Evaluation of global precipitation data sets over the Iberian Peninsula J. Geophys. Res. 2011 116 D20 10.1029/2010JD015481 94. Nashwan MS Shahid S Symmetrical uncertainty and random forest for the evaluation of gridded precipitation and temperature data Atmos. Res. 2019 230 104632 10.1016/j.atmosres.2019.104632 95. Waldman, D. Large-scale process-oriented modelling of soil erosion by water in complex watersheds (https://www.semanticscholar.org/paper/Large-scale-process-oriented-modelling-of-soil-by-Waldmann/19b903e68a26699d11eaaa314661e84473a28d19, 2010). 96. Diodato N Bellocchi G Storminess and Environmental Change 2014 Dordrecht Springer 97. Van Delden A The synoptic setting of thunderstorms in western Europe Atmos. Res. 2001 56 89 110 10.1016/S0169-8095(00)00092-2 98. Lim KJ Effects of calibration on L-THIA GIS runoff and pollutant estimation J. Environ. Manag. 2006 78 35 43 10.1016/j.jenvman.2005.03.014 99. Diodato N Filizola N Borrelli P Panagos P Bellocchi G The rise of climate-driven sediment discharge in the Amazonian River Basin Atmosphere 2020 11 208 10.3390/atmos11020208 100. Trouet Persistent positive North Atlantic Oscillation mode dominated the medieval climate anomaly Science 2009 324 78 80 10.1126/science.1166349 19342585 101. Guzzetti F Stark CP Salvati P Evaluation of flood and landslide risk to the population of Italy Environ. Manag. 2005 36 15 36 10.1007/s00267-003-0257-1 102. Gaume E A compilation of data on European flash floods J. Hydrol. 2009 367 70 78 10.1016/j.jhydrol.2008.12.028 103. Diodato N Historical evolution of slope instability in the Calore River Basin, Southern Italy Geomorphology 2017 282 74 84 10.1016/j.geomorph.2017.01.010 104. Jones, P. D. & Harris, I. C. Climatic Research Unit (CRU): Time-series (TS) datasets of variations in climate with variations in other phenomena v3. NCAS British Atmospheric Data Centre, http://catalogue.ceda.ac.uk/uuid/3f8944800cc48e1cbc29a5ee12d8542d (2008).