Nazzareno Diodato

Fredrik Charpentier Ljungqvist
email: fredrik.c.l@historia.su.se

Gianni Bellocchi

Fingerprint of climate change in precipitation aggressiveness across the central Mediterranean (Italian) area

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.

 .

 .

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 century [START_REF] Lionello | Objective climatology of cyclones in the Mediterranean region: A consensus view among methods with different system identification and tracking criteria[END_REF] and even earlier back to mid-seventeenth century on a restricted regional basis [START_REF] Camuffo | Improved understanding of past climatic variability from early daily European instrumental sources[END_REF] , but also from time-series of millennium-long hydrological extremes derived from documentary sources [START_REF] Diodato | A millennium-long reconstruction of damaging hydrological events across Italy[END_REF] . Long series derived from documentary data can help linking local and regional hydrological extremes to impacts [START_REF] Pichard | Hydroclimatology of the Lower Rhône Valley: historical flood reconstruction (AD 1300-2000) based on documentary and instrumental sources[END_REF] , assess climatic forcing features [START_REF] Benito | Recurring flood distribution patterns related to short-term Holocene climatic variability[END_REF] , and recognise the climatic variability in hazard-exposed areas [START_REF] Allan | Toward integrated historical climate research: The example of atmospheric circulation reconstructions over the earth[END_REF][START_REF] Diodato | Communicating hydrological hazard-prone areas in Italy with geospatial probability maps[END_REF] . 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 changes [START_REF] Glur | Frequent floods in the European Alps coincide with cooler periods of the past 2500 years[END_REF] .

In the Mediterranean region, local to regional climate experiences a mix of gradual and abrupt shifts of precipitating systems, including erratically distributed erosive storms [START_REF] Diodato | Decadal modelling of rainfall-runoff erosivity in the Euro-Mediterranean region using extreme precipitation indices[END_REF] . In fact, while northern Europe is primarily influenced by the North Atlantic storm tracks [START_REF] Rogers | North Atlantic storm track variability and its association to the North Atlantic oscillation and climate variability of Northern Europe[END_REF] , southern Europe is situated at the junction of major air-mass circulations, periodically arising from Atlantic, Mediterranean and Siberia influences [START_REF] Xoplaki | Modelling climate and societal resilience in the Eastern Mediterranean in the last millennium[END_REF] . 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 processes [START_REF] Dayan | Atmospheric conditions inducing extreme precipitation over the eastern and western Mediterranean[END_REF] . 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 rainfall [START_REF] Hofstätter | new classification scheme of European cyclone tracks with relevance to precipitation[END_REF] .

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 estimates [START_REF] Diodato | Decadal and century-long changes in the reconstruction of erosive rainfall anomalies at a Mediterranean fluvial basin[END_REF] , both in time (i.e., with annual resolution or finer) and space (i.e., non-locally calibrated) [START_REF] Poirier | Comparison of estuarine sediment record with modelled rates of sediment supply from a western European catchment since 1500[END_REF] . 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 westerlies [START_REF] Pinto | Extreme precipitation events over northwest Italy and their relationship with tropical-extratropical interactions over the Atlantic[END_REF] , fed by heat and moisture fluxes from the Mediterranean Sea and maintained by mesoscale processes which, in turn, determine a suite of distinct, spatiallyscattered erosive patterns (as in the schematic representation of Fig. 2).

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 scales [START_REF] Wainwright | Weathering, soils and slope processes[END_REF] . 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 records [START_REF] López-Vicente | Identifying erosive periods by using RUSLE factors in mountain fields of the Central Spanish Pyrenees[END_REF][START_REF] Borrelli | Rainfall erosivity in Italy: A national scale spatio-temporal assessment[END_REF][START_REF] Hernando | Estimate of the (R)USLE rainfall erosivity factor from monthly precipitation data in mainland Spain[END_REF] or non-uniform time intervals across different spatial scales [START_REF] Ballabio | Mapping monthly rainfall erosivity in Europe[END_REF] . 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 MedCA [START_REF] Diodato | Discovering historical rainfall erosivity with a parsimonious approach: A case study in Western Germany[END_REF] .

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. [START_REF] Diodato | A millennium-long reconstruction of damaging hydrological events across Italy[END_REF] . 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 R 2 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).

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 R 2 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 shift [START_REF] Wanner | North Atlantic oscillation-concepts and studies[END_REF] (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 statistics [START_REF] Alexandersson | Homogenization of Swedish temperature data. Part I: homogeneity test for linear trends[END_REF][START_REF] Pettitt | A non-parametric approach to the change-point problem[END_REF][START_REF] Buishand | Some methods for testing the homogeneity of rainfall records[END_REF][START_REF] Alexandersson | A homogeneity test applied to precipitation data[END_REF][START_REF] Worsley | Confidence regions for a change-point in a sequence of exponential family random variables[END_REF][START_REF] Wang | Penalized maximal t test for detecting undocumented mean change in climate data series[END_REF] 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 CE [START_REF] Miller | Abrupt onset of the Little Ice Age triggered by volcanism and sustained by sea-ice/ocean feedbacks[END_REF] ) 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) [START_REF] Wagner | The influence of volcanic, solar and CO 2 forcing on the temperatures in the Dalton minimum 1790-1830: A model study[END_REF] . 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 documented [START_REF] Trigo | The North Atlantic oscillation influence on Europe: Climate impacts and associated physical mechanisms[END_REF][START_REF] Pasini | Attribution of precipitation changes on a regional scale by neural network modeling: A case study[END_REF] . 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 patterns [START_REF] Hurrell | Decadal trends in the North Atlantic oscillation: Regional temperatures and precipitation[END_REF] . 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 variability [START_REF] Chelliah | Tropical multidecadal and interannual climate variability in the NCEP-NCAR reanalysis[END_REF] . 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 NAO [START_REF] Halpert | Climate assessment for[END_REF] . 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.

We derived that the interannual variability of rainfall erosivity is pronounced either before or after the changepoint (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 Mediterranean [START_REF] Angulo-Martínez | Do atmospheric teleconnection patterns influence rainfall erosivity? A study of NAO, MO and WeMO in NE Spain, 1955-2006[END_REF] ; and for Montenegro, on the east side of the MedCA 61 ). 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 erosivity [START_REF] Glur | Frequent floods in the European Alps coincide with cooler periods of the past 2500 years[END_REF] .

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 data [START_REF] Zolina | Improving estimates of heavy and extreme precipitation using daily records from European rain gauges[END_REF] . 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 subregional 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 season [START_REF] Diodato | How the aggressiveness of rainfalls in the Mediterranean lands is enhanced by climate change[END_REF] . 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 scales [START_REF] Diodato | Recalling the past of erosive rainfall hazard: A tandem with the future projections?[END_REF] . According to Cislaghi et al. [START_REF] Cislaghi | Statistical assessment of trends and oscillations in rainfall dynamics: Analysis of long daily Italian series[END_REF] and Pavan et al. [START_REF] Pavan | Daily precipitation observations over Emilia-Romagna: Mean values and extremes[END_REF] , 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 Tosti 67 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 century [START_REF] Myhre | Frequency of extreme precipitation increases extensively with event rareness under global warming[END_REF][START_REF] Rädler | Frequency of severe thunderstorms across Europe expected to increase in the 21st century due to rising instability[END_REF] .

Climate model simulations with the new generation of Coordinated Downscaling Experiment over Europe (EURO-CORDEX [START_REF] Jacob | EURO-CORDEX: New high-resolution climate change projections for European impact research[END_REF] ) 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 precipitation [START_REF] Alfieri | Global warming increases the frequency of river floods in Europe[END_REF] . 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. [START_REF] Alpert | The paradoxical increase of Mediterranean extreme daily rainfall in spite of decrease in total values[END_REF] , Paxian et al. [START_REF] Paxian | Present-day and future mediterranean precipitation extremes assessed by different statistical approaches[END_REF] and Caloiero et al. [START_REF] Caloiero | The trend of monthly temperature and daily extreme temperature during 1951-2012 in New Zealand[END_REF] but questioned by Mariani and Parisi [START_REF] Mariani | Extreme rainfalls in the Mediterranean area. In: Storminess and environmental change-Climate forcing and responses in the Mediterranean region[END_REF] for the whole basin and the western sub-basin. Other researchers 75 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 rainstorms [START_REF] Panagos | Spatio-temporal analysis of rainfall erosivity and erosivity density in Greece[END_REF] .

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 climate [START_REF] Liu | Spatial-temporal changes of rainfall erosivity in the loess plateau, China: Changing patterns, causes and implications[END_REF] . 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 states [START_REF] Zhang | Complexity in estimating past and future extreme short-duration rainfall[END_REF] . 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:

1. 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. 2. 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. 3. The way how in the recent decades Mediterranean cyclones have been producing trends of erosivity from rainfall extremes remains elusive [START_REF] Lionello | Objective climatology of cyclones in the Mediterranean region: A consensus view among methods with different system identification and tracking criteria[END_REF][START_REF] Camuffo | Improved understanding of past climatic variability from early daily European instrumental sources[END_REF] . 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. [START_REF] Berg | Strong increase in convective precipitation in response to higher temperatures[END_REF] 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 warming [START_REF] Aghakouchak | How do natural hazards cascade to cause disasters?[END_REF] , though with distinct patterns of change in small domains without the emergence of spatially and temporally homogeneous trends [START_REF] Libertino | Evidence for increasing rainfall extremes remains elusive at large spatial scales: The case of Italy[END_REF] .

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 reliefs [START_REF] Diodato | a rainfall erosivity model for the Mediterranean region[END_REF] . It was also shown [START_REF] Diodato | Assessing and modelling changes in rainfall erosivity at different climate scales[END_REF] 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 rates [START_REF] Diodato | Drought stress patterns in Italy using agro-climatic indicators[END_REF] . The regional analysis may also obscure sub-regional trends [START_REF] Zittis | Observed rainfall trends and precipitation uncertainty in the vicinity of the Mediterranean, Middle East and North Africa[END_REF] . 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) [START_REF] Diodato | Decadal and century-long changes in the reconstruction of erosive rainfall anomalies at a Mediterranean fluvial basin[END_REF] and in the Po river basin (northern Italy) [START_REF] Diodato | Historical predictability of rainfall erosivity -a new reconstruction for monitoring extremes over Northern Italy (1500-2019 CE)[END_REF] , showed similar increasing trends since the end of the LIA. Diodato et al. [START_REF] Diodato | 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[END_REF] 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 region [START_REF] Borrelli | An assessment of the global impact of 21st century land use change on soil erosion[END_REF] . 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 stormrelated 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 Diodato [START_REF] Diodato | Predicting RUSLE (Revised Universal Soil Loss Equation) monthly erosivity index from readily available rainfall data in Mediterranean area[END_REF] and Acquaotta et al. [START_REF] Acquaotta | Estimation of rainfall erosivity in Piedmont (Northwestern Italy) by using 10-minute fixed-interval rainfall data[END_REF] , with updating from the network of SCIA (http://www.scia.ispra mbien te.it/wwwro otsci a/scia.html)-National System for the Collection and Elaboration of Climatological Data [START_REF] Desiato | Climate indicators for Italy: Calculation and dissemination[END_REF] , 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 ://tinyu rl.com/h4juz uv) 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).

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 resolution [START_REF] Pauling | Five hundred years of gridded high-resolution precipitation reconstructions over Europe and the connection to large-scale circulation[END_REF] . The GPCC (Global Precipitation Climatology Centre) provides an extension of the seasonal precipitation dataset until 2019 (retrieved from https ://tinyu rl.com/wlnvz nq). GPCC products are known to outperform other similar products (e.g. CRU, ERA-Interim, ERA-40) in Mediterranean areas [START_REF] Belo-Pereira | Evaluation of global precipitation data sets over the Iberian Peninsula[END_REF][START_REF] Nashwan | Symmetrical uncertainty and random forest for the evaluation of gridded precipitation and temperature data[END_REF] .

A categorical variable, the annual severity storm index sum (SSI) from Diodato et al. [START_REF] Diodato | A millennium-long reconstruction of damaging hydrological events across Italy[END_REF] 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 approach [START_REF] Renard | Predicting soil erosion by water: a guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE)[END_REF] . The non-linear model takes the following form: 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 ( P Sum • P Aut ) 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 ( P Sum • P Aut ) was considered here because the power of rainfall exerts the maximum of erosive forces in summer and autumn (Fig. 6).

'0"E 18°0'0"E 15°0'0"E 15°0'0"E 12°0'0"E 12°0'0"E 9°0'0"E 9°0'0"E 45°0'0"N 45°0'0"N 42°0'0"N 42°0'0"N 39°0'0"N 39°0'0"N 1 2 3
Severe convective erosive storms occur between June and November in Europe [START_REF] Diodato | Decadal modelling of rainfall-runoff erosivity in the Euro-Mediterranean region using extreme precipitation indices[END_REF] . 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 frequent [START_REF] Diodato | Decadal modelling of rainfall-runoff erosivity in the Euro-Mediterranean region using extreme precipitation indices[END_REF] . 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 (P Win + P Spr ) 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 totals [START_REF] Diodato | A millennium-long reconstruction of damaging hydrological events across Italy[END_REF] . The additive component P Win + P Spr supports the assumption of a softer dependence of erosivity on precipitation intensity.

The concept of the model is summarised in Fig. 7 following Waldam [START_REF] Waldman | Large-scale process-oriented modelling of soil erosion by water in complex watersheds[END_REF] and Diodato and Bellocchi [START_REF] Diodato | Storminess and Environmental Change[END_REF] . 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 flooding [START_REF] Van Delden | The synoptic setting of thunderstorms in western Europe[END_REF] . 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 SD Max(Ps) . In this case, raindrops with high kinetic energy may cause local downpours.

Trend assessment of extreme values. The quantile approach was applied to identify step trends at any ED X 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:

Figure 1 .

 1 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.europ a.eu/ conte nt/globa l-rainf all-erosi vity) 55 .

Figure 2 .

 2 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).

Figure 3 .

 3 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). (a 1 ) 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).

Figure 4 .

 4 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/getpa ge.pl).

  10:22062 | https://doi.org/10.1038/s41598-020-78857-3 www.nature.com/scientificreports/

Figure 5 .

 5 Figure 5. Spatial distribution of digital rain gauges used for calibration (black numbers) and validation (red numbers) across the Mediterranean Central Area.

Figure 6 .

 6 Figure 6. Seasonal regime of storm cells with the related scheme of variability of the monthly storm erosivity strength over the MedCA.

( 1 )

 1 REMM = A • SD Max(Ps) • P Sum • P Aut + B • α + ASSIS(GF) k • P Win + P Spr + C (5) Q(ED Xt-w ) T = exp µ ED * Xt-w + u T • σ ED * Xt-w

Figure 7 .

 7 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 Dataset 104 . The terms of Eq. (1) are reported.

Scientific Reports | (2020) 10:22062 | https://doi.org/10.1038/s41598-020-78857-3

Scientific Reports| (2020) 10:22062 | https://doi.org/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.

Funding

Open Access funding provided by Stockholm University.

where Q(ED Xt-w ) T is the j th rainfall erosivity density (ED X ) quantile of the log-normal distribution with assigned return period; µ ED * Xt-w and σ ED * Xt-w are the mean and the standard deviation of the variable ED X t * = ln(ED X t ). The subscript t-w indicates the computation of the generic variable at the t-time over a 22-year moving window (w); u T 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 changes [START_REF] Alfieri | Global warming increases the frequency of river floods in Europe[END_REF] .

Model calibration and assessment.

To assess the model, statistical analyses were performed with STAT GRA PHICS (http://www.duke.edu/~rnau/sgwin 5.pdf), with the graphical support of WESSA (https :// www.wessa .net) and CurveExpert routines (https ://www.curve exper t.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 ≤ R 2 ≤ 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 uncertainty [START_REF] Lim | Effects of calibration on L-THIA GIS runoff and pollutant estimation[END_REF] . 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 R 2 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 achieved [START_REF] Diodato | The rise of climate-driven sediment discharge in the Amazonian River Basin[END_REF] . 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.

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.europ a.eu/conte nt/globa l-rainf all-erosi vity. 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/getpa ge.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://clime xp.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.ispra mbien te.it/wwwro otsci a/scia.html, completed with the data provided from the digital National Agrometeorological Network accessible through https ://tinyu rl.com/h4juz uv. Updated seasonal precipitation datasets were retrieved from the Global Precipitation Climatology Centre (GPCC) through https ://tinyu rl.com/wlnvz nq. 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 stormseverity index inputs), are available in the supplementary Table S1.

Author contributions

Competing interests

The authors declare no competing interests.