Arno is the second most important river of central Italy (241 km long). It originates on Mount Falterona at 1385 m a.s.l. and flows into the Tyrrhenian Sea at Marina di Pisa (43°40′N, 10°16′E). The river flooded Florence quite often in historical times, most recently in 1966 (November 4), with 4500 m³ s⁻¹ after rainfall of 437.2 mm in Badia Agnano (43°26′N, 11°38′E, 271 m a.s.l.) and 190 mm in Florence, in 24 h (Martin, 1967). Located across the northern Apennine chain (Figure 3), the Arno River Basin (ARB), from its origin to the outlet, has an area of ∼8229 km², with an average elevation of ∼353 m a.s.l. The Casentino forest massif, in the province of Arezzo (43°28′N, 11°52′E, 134–954 m a.s.l.), is the valley in which the first tract of the Arno River flows. Then, it enters the upper Valdarno (highlighted in Figure 3), a long valley bordered on the east by the Pratomagno massif and on the west by the hills around Siena (43°19′N, 11°19′E, 200–510 m a.s.l.).

In the middle and lower Valdarno, the course of the Arno forms an important natural demarcation between the areas north of the river, which are cooler and rainier, and those south, which are warmer and more arid. The climate of the ARB is determined by its position in relation to the direction of wind motion (with an intense western disturbance) and geographical factors such as distance from the sea, and elevation, orientation and distribution of the reliefs, which determine a marked differentiation of climates. The highest mean annual precipitation values, >1500 mm, occur on the crests of Pratomagno and the upper Casentino. The Arno receives water from several tributaries. From the right side, the river Sieve comes from Mugello reliefs (43°95′N, 11°38′E, 162–437 m a.s.l.), while in the southern part of the basin (Val di Chiana, 43°15′N, 11°49′E, 405 m a.s.l. on average), the rivers Pesa and Era form the valleys of the left-side sub-basins. Here, annual precipitation varies between 800 and 900 mm. Rainfall intensity >100 mm d⁻¹ occurs in most of the basin, with peaks of 180–200 mm d⁻¹ in some sectors of the Apennines (Monti and Rapetti, 2011).

Annual (R)USLE-based erosivity data were obtained for 1964–1990. This calibration period covers two climatic regimes, with a temperature decrease and weak winter westerlies until about 1976, and a strong temperature increase and winter westerlies after about 1976 (Mariani and Diodato, 2014). The actual erosivity data were obtained from Bazzoffi and Pellegrini (1990) at Vicarello di Volterra and Cevasco et al. (2015) at Levanto, and derived from Borselli et al. (2004) and Angeli et al. (2007) at Florence based on a power relationship between the R-factor and the monthly rainfall data (source: SCIA—National System for Climate Data Collection and Dissemination, http://www.scia.isprambiente.it). With about 2300 (Levanto), 1718 (Florence) and 731 (Vicarello di Volterra) MJ mm hm⁻² h⁻¹ yr⁻¹ on average, the three sites (Figure 1C, red circles) are representative of situations occurring in three sectors of the ARB. Located on the Ligurian coast (north-western Mediterranean) and at the foothills of Tuscany, Levanto is representative of hydrological conditions occurring in upper Valdarno, characterized by heavy precipitation events leading to large amounts of rainfall (Hofstätter et al., 2016). Florence is located in the southwest corner of the middle Valdarno, marked by annual precipitation totals of 800–900 mm. Vicarello di Volterra is located downstream of Florence, in the lower Valdarno, a broad valley where the average annual rainfall is below 800 mm. The weighted areal mean of actual rainfall erosivity (MJ mm hm⁻² h⁻¹ yr⁻¹) for the Arno Basin (RS) was calculated for each year of the period 1964–1990 as a linear combination of the site-specific erosivity values (S₁: Levanto; S₂: Florence; S₃: Vicarello di Valdarno): R S = s ⋅ ( a 1 ⋅ S 1 +   a 2 ⋅ S 2 +   a 3 ⋅ S 3 ) (1) where one scale factor s = 1.35 and three relative weighting factors a ₁ = 0.2, a ₂ = 0.1 and a ₃ = 0.7 are calibrated values to converge to the long-term areal mean (1541 MJ mm hm⁻² h⁻¹ yr⁻¹; Figure 1B). The estimated scale factor 1 < s < 2 is an upscaling parameter, which is consistent with the assumption that individual sites are fractal elements of the area (e.g., Diodato et al., 2017b). The estimated relative weighting factors a _i (i = 1, 2, 3) reflect the characteristics of the three stations (e.g., Bartolini et al., 2018). At the low end of the ARB, Vicarello di Volterra (station S₃) has the greatest influence on the areal erosivity because it is located upwind of the western currents that are the bearers of heavy rainfall. As they advance eastward and northeastward, many of these storms transfer energy in the form of erosive rain over much of the basin. Outside the basin, Levanto (station S₁) also maintains an upwind position with respect to the northern Apennine slopes, but thunderstorms are forced to move east and northeast, thus passing only over some areas of the river basin. Finally, Florence (station S₂) contributes the least as it is located far from the sea and is not much affected, on average, from the heaviest rainfall that discharges its energy before reaching the easternmost part of the basin.

We used Eq. 2 to reconstruct a time-series of annual storm erosivity (MJ mm hm⁻² h⁻¹ yr⁻¹) in the ARB from 1000 to 2019 CE. We also obtained erosivity density—Eq. 3—using the rainfall (mm) data provided for the period 1891–2016 (∼900 mm on average) by the GPCC (Global Precipitation Climatology Center) by 0.25° monthly land-surface precipitation dataset Schneider et al. (2018), built on global telecommunication system-based and historical data (https://opendata.dwd.de/climate_environment/GPCC/html/fulldata-monthly_v2018_doi_download.html).

In addition to extensive research in the main Italian historical archives and libraries, literary sources were consulted by web search (https://books.google.com), which generated ∼100,000 bibliographical documents. We collected a massive bibliographic dataset but only ∼500 records have met the criteria of including the keywords abundant rainfall, storm, downpour, diluvial, flood and alluvial (piogge abbondanti, tempesta, nubifragio, diluvio, inondazione, alluvione), as well as some Latin locutions (e.g., magnae pluviae, aqua maxima, diluvium, excrescentia fluminum, inundatio) that were chosen for careful reading. Useful data were extracted from documentary sources by converting the information contained in historical accounts (Morozzi, 1766; Aiazzi, 1845; Giorgini, 1854; Rossini, 1855; Becchi and Paris, 1989; Salvestrini, 2005; Menduni, 2017; Pelucani et al., 2017; Ricci et al., 2017; Rombai et al., 2017; Salvestrini et al., 2017) into numerical values on an index scale (see Supplementary Material S1—erosivity reconstruction in Supplementary Material). Transformation process required a dynamic understanding of the historical information used in the analysis, with an in-depth knowledge of regional and sub-regional climates, and familiarity with the relative strengths and weaknesses of each type of source. Weather hind-casting (Pfister et al., 2018) was used to familiarize with well-documented anomalies in the instrumental period before analyzing similar cases in the pre-instrumental time. Storms occurring during the summer season (May to August) were not considered in this study, as they generally only affect small or isolated areas of the ARB and are not representative for long-term storm reconstruction. A scoring system was thus established (Wetter et al., 2011) to rate the MSSI (defined in the time of year between September and April) as 0 (normal), 1 (stormy), 2 (stormy with some flooding), 3 (stormy with important flooding), and 4 (extraordinarily stormy with major flooding):

Normal: average storm or storm passed unobserved, without comments about its severity or its impacts on society and economy.

The study was based on the systematic and critical analysis of data on the above phenomena provided by Italian documentary sources. For most of the information, it was possible to carry out an “event check” considering more than one documentary source on the same event. It was also possible to contextualize the storms with other types of historical events (e.g., social, agricultural, religious). In this way, the reliability of information was assessed going beyond quantitative data and looking for other sources of information such as diaries, newspapers, chronicles and local stories.

For the historical reconstruction of annual rainfall erosivity (MJ mm hm⁻² h⁻¹ yr⁻¹) in the ARB, we developed the Rainfall Erosivity Model for the Arno River Basin (REM _ARB), which combines low-resolution explanatory inputs to form a parsimonious concept compatible with the scenario depicted by (R)USLE-based erosivity data (Renard et al., 1997). The non-linear model of annual erosivity in the ARB takes the following form: R E M ARB = A ⋅ [ α + ( ∑ i ∈ X MSSI i ) k ] + B ⋅ [ β + S D i ∈ Y ( MSSI i ) k ] ⋅ ( γ + ∑ i ∈ Z MSSI i ) (2) where: A (MJ mm⁻¹ hm⁻² h⁻¹ yr⁻¹) and B (MJ mm⁻¹ hm⁻² h⁻¹ yr⁻¹) are scale parameters converting the result of two dimensionless terms of the model into the output unit; α, β and γ are shift parameters estimating rainfall erosivity when the monthly input values (i = 1, …, 12 months) of the storm-severity index (MSSI_i) are all equal to zero; k is a shape parameter; SD is standard deviation; X (i = 9, 10, 11), Y (i = 1, 2, 3, 4, 9, 10, 11, 12) and Z (i = 1, 2, 3, 4, 12) are different groups of months (where the symbol ∈   stands for “belongs to”). The erosivity model was designed not to contain a position parameter, i.e., a constant term prescribing the fixed amount that contributes to rainfall erosivity when all other terms of the equation are equal to zero. In this way, the model predicts a null rainfall erosivity for null input values.

The concept of the model is summarized in Figure 4. Complex climatic characteristics affect erosive rainfall, depending on rain-splash, runoff and flood generation mechanisms, and it is difficult to disentangle the climatic component from substantial natural variability and anthropogenic impacts (Sofia and Nikolopoulos, 2020). According to Waldman (2010) and Diodato and Bellocchi (2014), non-linear relationships emerging across river basins depend on the processes that dominate a given hydrological regime.

In southern Europe, erosive storms mostly occur in autumn, when the occurrence of severe thunderstorms is intimately associated with convective processes (e.g., Guzzetti et al., 2003; Diodato and Bellocchi, 2012; Petrucci et al., 2014; Ballabio et al., 2017). Spring precipitation likely leads to a different regime of hydrological extremes compared to September–October, when extreme floods are more frequent. This indicates that the period September–April delineates a key time-window in the year for estimating rainfall erosivity in fluvial basins. The separation between advective and convective events is considered important, assuming that convective precipitation (short and intense) has a greater dynamic and variability than generally more static advective events.

The seasonal differences in precipitation are governed by convective processes, more dynamic and frequent in summer (when rain-splashes dominate), or by a variable mixture of convective and advective precipitation, more frequent in autumn (when overland flows dominate). Then, winter and spring rainfall (long-lasting and low-intensity), which usually originates from broad frontal activity, carries large volumes of rainwater (through orographic stratiform precipitation) that cause large-scale hydrological processes like flooding (Van Delden, 2001). These processes are mostly captured in the REM _ARB—Eq. 2—by MSSI. Very intense convective events, or those of a mixed advective-convective nature (with high kinetic energy causing local showers and leading to splash-erosivity processes) are interpreted by the SD term. The erosivity density (ED, MJ hm⁻² h⁻¹) in the ARB (ED _ARB) was calculated for the 1891–2016 period as: E D ARB = R E M ARB P ARB (3) where P _ARB (mm yr⁻¹) is the annual precipitation in the basin.

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 AgriMetSoft Online Calculators (https://agrimetsoft.com/calculatorshttps://www.curveexpert.net/). The parameters of Eq. 2 were calibrated against actual rainfall erosivity data according to statistical criteria. The first condition was to minimize the distance between modeled and actual erosivity data, by minimizing the Mean Absolute Error (optimum, 0 ≤ MAE < ∞, MJ mm hm⁻² h⁻¹ yr⁻¹). 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 maximize the determination coefficient (0 ≤ R² ≤ 1, optimum) that is the variance explained by the model. The third conditions approximates the unit slope of the straight line that would minimize the bias of the linear regression actual vs. modeled 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 index (−∞<EF ≤ 1, optimum; Nash and Sutcliffe, 1970) has also been calculated as an indicator of model performance uncertainty because values above 0.6 indicate limited model uncertainty, presumably associated with narrow parameter uncertainty (Lim et al., 2006). 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 modeling solutions with small MAE and large R ² values were obtained. Then, for the final selection, a third criterion −|b−1| = min − was additionally involved. Each predictor was repositioned over >50 iterations until convergence was achieved. 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 observation-based data.

We have extracted 369 extreme hydrological events, occurring in the ARB from 1000 to 2019 CE. The classification of these events by severity led to 213 stormy events, 98 very stormy events, 41 great stormy events, and 17 extraordinary stormy events. However, our historical hydrological database is based on several types of heterogeneous sources, including reports that may be an exaggeration of what was an uncritical reference to previous sources, printing errors in the documentation, and natural records of environmental factors (“proxy” data). All of them contain different sources of uncertainty in the “cataloguing” of storm records (Pavese et al., 1994). It is well ascertained, for instance, that small storms tend to be underestimated, especially those far from inhabited centers, while more localized ones can be frequent. In order to resolve some of these uncertainties in our database, we have established a reasonable standard for Monthly Storm-Severity Index (MSSI) recorded events. This was facilitated by defining a partition of the time-series in three equally-long sub-periods following Diodato et al. (2019b),—1000–1340 (Figure 5A), 1341–1680 (Figure 5B), and 1681–2019 (Figure 5C)—and verifying, for each sub-period and for the entire dataset (Figure 5D), the scale-invariance in the relationship between the number of events larger than the strength of a given storm event and the same strength event. The completeness analysis was formalized with the relationship between the cumulative number of events (CEN) and the MSSI values in the range 1 ≤ MSSI ≤4, as follows: log 10 ( C E N i , j ) = a + b ⋅ MSSI i , j     with   i = 1 , … ,   4     and   j = 1 , … ,   4 (4) where MSSI is the Monthly Storm-Severity Index by severity class (i) and sub-period (j). In total, 369 events within the range 1 ≤ MSSI ≤4 are described in qualitative terms as stormy, very stormy, great stormy and extraordinary stormy.

The negative slopes in the three sub-periods, and for the entire dataset, reflect the principle of a progression toward a lower frequency as storms become more severe. With determination coefficients 0.97 ≤ R ² ≤ 1.00, it can be assumed that storms in the period 1000–2019 CE are scale-invariant and that ancient sources from different periods are not too fragmentarily available to facilitate a quantitative reconstruction (Mazzarella and Diodato, 2002).

To estimate the annual values of the mean rainfall erosivity in the ARB, we developed and calibrated a simplified statistical model (REM _ARB)—Eq. 2—which summarizes (through MSSI inputs) the relationship between spatial patterns of climate and storm erosivity, consistent with a sample (1964–1990) of detailed (R)USLE-based data obtained for the study area with Eq. 1. The mean actual erosivity is 1541 MJ mm hm⁻² h⁻¹ yr⁻¹ (±340 standard deviation), the values varying between 1134 (in 1977) and 2306 (in 1981) MJ mm hm⁻² h⁻¹ yr⁻¹. In the same period (1964–1990), the calibrated model gave 1542 MJ mm hm⁻² h⁻¹ yr⁻¹ on average (±340 standard deviation), and 1256 (in 1977) and 2487 (in 1966) MJ mm hm⁻² h⁻¹ yr⁻¹ as minimum and maximum estimates, respectively. A qualitative validation for the 1951–2010 period was also introduced in order to control if modeled rainfall erosivity can be representative of the erosional sediment processes that happen across the ARB. For the calibration period 1964–1990, we obtained the coefficients A = 58.85 MJ mm⁻¹ hm⁻² h⁻¹ yr⁻¹, B = 1.00 MJ mm⁻¹ hm⁻² h⁻¹ yr⁻¹, α = 10.00, β = 4.61, γ = 2.46 and k = 2.00 in Eq. 2. With these values, the linear regression between actual and estimated erosivity data is statistically significant (F test p∼0.00). The R ² statistic (goodness of fit) indicates that the REM _ARB explains 88% of the erosivity variability. MAE (mean absolute error) equal to 96 MJ mm hm⁻² h⁻¹ yr⁻¹ and MAPE (mean absolute percent error) equal to 7%, with Kling-Gupta Efficiency (KGE) equal to 0.91, indicate satisfactory model performance and efficiency. Although there is an indication of possible serial correlation in the residuals (Durbin-Watson statistic p < 0.05), the calibrated regression (Figure 6A, black dotted line) shows only negligible departures of data-points from the 1:1 identity line (red line), which indicates the substantial ability of the model, Eq. 2, to predict actual erosivity. In fact, the Nash-Sutcliffe efficiency value (EF = 0.88) indicates limited uncertainty in model estimates. In particular, the regression line has an intercept a = −0.58 (±114.4) MJ mm hm⁻² h⁻¹ yr⁻¹ and a slope b = 1.00 (±0.1) near or equal to the optimum values (a = 0 and b = 1). Figure 6B indicates the normal approximation of model residuals (normality test p > 0.05, Jarque and Bera, 1981). The distribution of the percentiles of modeled rainfall erosivity (Figure 6C, orange curve) approaches the distributional shape of the observed erosivity, indicating a satisfactory prediction over the range of erosivity values (including the low and high erosivity values).

Di poi le molte pioggie, accrescimento dei fiumi, con ispessi lavamenti à dispogliati in parte l’alte cime d’essi monti, lasciando il loco della terra, il sasso si trova essere circondato dall’aria, e la terra d’essi lochi partita

The qualitative validation shows that modeled erosivity can be a representative driver of erosional sediment processes occurring in the ARB (Figure 6D). The interannual variation in the volumes of soil potentially transported through the ARB shows an overall downward trend, although with some large fluctuations, and this is sufficiently described by our modeled erosivity. Some discrepancies between the two lines could be due to inaccurate determinations of the transport of sand in the Arno river bed, in particular for the year 1960, for which no extreme rainfall was detected in any of the ARB pluviometric station. According to Monti and Rapetti (2011), the extraction of large volumes of inert materials can be associated with an intense potential transport, concentrated between the late 1950s and the first half of the 1960s. We have confirmation that pluviometric events (not soil and landscape features) are the main forcing agent of sediment from authors’ own words (Monti and Rapetti, 2011, p. 124):

This analysis then points out that the baseline categorical MSSI was a reliable substitute of storm rainfall.

Figure 7A shows the evolution of the estimated annual rainfall erosivity values (blue circles) obtained by Eq. 2 for the ARB (areal mean) over the 1,000–2019 CE period: mean 1,412 (±393 standard deviation), minmum 1256, maximum 4892 MJ mm hm⁻² h⁻¹ yr⁻¹. The distributional properties of the actual data sample and the millennium-long reconstructed time-series are more similar as moving toward the highest percentiles (data not shown). An issue is how to develop a continuous time-series of the annual output when integer monthly inputs are present. In our case, the use of integer values for the MSSI may not sufficienly reflect the fluctuations that occur in the annual rainfall erosivity. For a better identification of possible trends and oscillations in the estimated data and compare contemporary and historical patterns, the time-series of annual rainfall erosivity was also filtered using a 11-years low-pass Gaussian function (Figure 7A, bold blue curve). We observe that the world’s longest series of annual erosive rainfall, dating back to the latter half of the Medieval Climate Anomaly (MCA; here 1000–1249 CE), is associated with frequent low-erosivity storms (∼1300 MJ mm hm⁻² h⁻¹ yr⁻¹), which occurred once every four years. During the MCA, infrequent extreme erosivity values (∼2000 MJ mm hm⁻² h⁻¹ yr⁻¹) were only present in ∼3% of the years (i.e., in 1011, 1029, 1050, 1060, 1080, 1086, 1117, 1050, 1129, 1139, 1272, 1282 and 1333), while the years 1168 and 1179 are exceptional, with ∼4000 MJ mm hm⁻² h⁻¹ yr⁻¹. For these years, the 12th-century chronicler Bernardo Maragone reported nine and 13 floods, respectively (Salvestrini et al., 2017). This long-term erosivity pattern is similar to that found by Longman et al. (2019) for the eastern Mediterranean region, reflecting less rainy conditions with limited hydroclimatic variability during much of the MCA. With the onset of the Little Ice Age (LIA; here 1,250–1849 CE), erosivity often became more vivid (likely associated with the enhanced variability of dynamical processes, e.g., Raible et al., 2018), with values mostly just below 2000 MJ mm hm⁻² h⁻¹ yr⁻¹ (corresponding to the 95th percentile of the whole series, yellow line in Figure 7A). These natural events had disastrous consequences in much of western Europe (Pouzet and Maanan, 2020). The second decade of the 14th century was affected by exceptionally high rainfall in Atlantic Europe, which had “water at ankles” level, according to the iconic formula of Le Roy Ladurie (1982). Taking special care in describing damaging hydrological phenomena over the 14th century, the chroniclers’ texts convey useful information not only to know the climate of that century but also the characteristics of climate impacts on society, agriculture and transport, and the measures taken to protect the territory surrounding the river and the town of Florence. Giovanni Villani (c. 1280–1348) reports the catastrophic flooding associated with the erosive event of 1333 (estimated value equal to 2509 MJ mm hm⁻² h⁻¹ yr⁻¹), whose chronicle is an example of the importance of the narrative form used in historical accounts (as from the posthumous work Villani, 1587, p. 673): ... per lo modo che chi lo leggerà per lo tempo avvenire, potrà comprendere i termini fermi e notabili onde faremo menzione appresso (“... for the way that those who will read it for the time to come, will be able to understand the firm and notable terms that we will mention below”).

That exceptional event was the result of a more general regional-scale shift, which can be considered the main cause of the hydrological change that occurred during that period, also in other areas of northern Italy (Glur et al., 2013; Baldini and Bedeschi, 2018; Diodato et al., 2020c). Leonardo da Vinci himself, who knew the territory of Tuscany well, left many drawings with projects that were supposed to alleviate the effects of the increasingly dangerous floods and soil erosion (Pelucani.et al., 2017). On sheet 126 υ. of the Codex Atlanticus Leonardo writes (De Lorenzo, 1871, p. 122):

A significant change-point in the erosivity time-series was detected in 1490 (red arrow in Figure 7A), in the central part of the LIA, either with the T-Test (Panofsky and Brier, 1958), the Worsley’s likelihood ratio (Worsley, 1986) or the Standard Normal Homogeneity Test (Alexandersson, 1986). After this change-point, erosivity has become more changeable, with erosive storms tending to oscillate more, as also shown by the 11-years Gaussian filtered curve (Figure 7A, bold blue curve). The statistically relevant change-point provides a picture of variations in erosivity related to changes in climate patterns, with cold conditions dominating during the LIA. Hence, we note groups of annual values around and above the 98th percentile in this cold period (blue circles near the horizontal bold pink line in Figure 7A).

The uncertainty of hydrological information poses challenges for the analysis of historical storm data because the aggressiveness of storms evolves in complex hydrological processes, with local interactions between environmental context, climate fluctuation and human impact (Bintliff, 2002). The driving forces and actors in a close and reciprocal interaction result in cause-and-effect relationships that include feedbacks, where the erosive hazard becomes an expression of the human influence and its impact on landscape change. For instance, there is a slope effect on radiative balance and convective events (e.g., Bois et al., 2020). Thus, storm records after c. 1490 may have become more changeable due to climatic variations and human influence on slopes, resulting in less water infiltration by continuous deforestation and consequently more water available for runoff, reaching in an easier way the water lines and leading to a reduction in the time of water concentration in the basin. After the change-point, the construction of canals and ports along the Arno river led to the continuous deforestation of the timber for the consequent construction of ships, houses and other buildings (Salvestrini et al., 2010). We know that mount Morello, near Florence, already appeared deforested in the 14th century (Hofman and Perulli, 2000). Then, during the 15th and 16th centuries, the progressive deforestation of the Apennine chain and the inland reliefs accentuated the risk of seasonal runoff and flooding (Salvestrini.et al., 2017). Also in the 19th century, Targioni-Tozzetti (1852) lamented the absence of forest containment in the Grand Duchy of Tuscany for many years. Thus, natural factors related to the torrential nature of the rivers (based on the concentration of rainfall during the winter period from September to April), the steep slope of the elevation profile and the long and continuous deforestation, as well as the climatic surge triggered by solar minima, must be attributed to the clearly erosive character of the storms in the ARB until the end of the Litle Ice Age. In fact, we know that between the end of the 15th century and throughout the 16th century, the Arno River Basin was characterized by an increase in storms and floods, perhaps even favored by a worsening of the climate with the resumption of equinoctial autumn storms. Among others, extraordinary phenomena occurred in the years 1494, 1532, 1543, 1547, 1557, 1571 and 1589, characterised by a particularly strong erosive forcing of estimated 3483, 2431, 2487, 2872, 2488, 2431 and 3391 MJ mm hm⁻² h⁻¹ yr⁻¹, respectively. These events were particularly dangerous at the beginning of autumn, with concomitant tilling of soils. However, despite the awareness of the causes of the damage (deforestation, steep slopes, exploitation of the river bed) or, more simply, the fear that disasters would happen again, the Florentine administrators did not have the power to implement any protection measure because of the fragmentation of authorities and competences within the Arno River Basin: the discourse was dominated by economic considerations (timber supply, food production, protection of factories within the walls) that overlooked environmental and safety concerns (Salvestrini et al., 2010). It is precisely during this period that the Arno River Basin became most vulnerable, as it was also subject to cyclically recurrent deforestation (Ricci et al., 2017; Berti et al., 2019). It is hardly a coincidence that the climate in the central and final part of the LIA provided some injuring scenarios, with storms of unprecedented strength. In particular, the most striking group of erosive events occurred between the late 17th century and late 18th century, when we estimated erosivity values of ∼1600 MJ mm hm⁻² h⁻¹ yr⁻¹ once every about four years, with estremes of ∼3000 MJ mm hm⁻² h⁻¹ yr⁻¹ once every ∼20 years. The two maximum estimated erosivity values for the whole series, 4892 and 3844 MJ mm hm⁻² h⁻¹ yr⁻¹, were reached in 1688 (out-of-scale value of Figure 7A) and 1687, respectively. These annual rates are considerable, considering that they are smoothed over the basin area, which means that even higher erosive forcing may have occurred locally in the basin. This phase has also recently been identified as erosive-prone, characterized by high runoff in the southern Carpathians (Longman et al., 2019) and flooding in most of Europe (Blöschl et al., 2020). It corresponds to the grand solar (Maunder) minimum of 1645–1715 CE (Eddy, 1976), during which sunspots became exceedingly rare with lower-than-average temperatures in Italy (Figure 5, bottom colored band). Maximum erosivity values occurred just with such low temperatures, with likely a greater thermal contrast between the cold fronts and the warmer Mediterranean waters, which exacerbated the thermal difference between the air and the sea, making the precipitation events more intense. An increase in high-magnitude floods was observed in central and southern Europe c. 1700, linked to the cold and dry climate of the late Maunder minimum (Mudelsee et al., 2004), a key segment of the LIA for studying decadal-scale climatic change in Europe (Luterbacher et al., 2001).

During the 20th century, erosive fluctuations began a new, albeit modest rise, with some notable values in 1938, 1966 and 2013, comparable to those of the LIA. By the onset of the 20th century, the agricultural landscape already experienced the general phenomenon of rural abandonment that become typical during the second half of the century, with a recovery of wooded areas and an urbanization of parts of the countryside to the detriment of agricultural and pasture lands (Berti.et al., 2019). Then, from 1952 onwards, the application of the mountain law led to a 15-year period of renewed fervor not only with reforestation but also with the construction of infrastructures, with the participation of peripheral bodies like municipal consortia and mountain communities (Hofman and Perulli, 2000). Summary statistics for three climatic periods (Table 1) show that although the mean rainfall erosivity has continued to increase from medieval times until the recent period of warming (1342–1441 MJ mm ha⁻¹ h⁻¹ yr⁻¹), it is during the LIA that the coefficient of variation (30%) is maximized and the 95th and 98th percentiles are particularly high (2,431 and 2760 MJ mm ha⁻¹ h⁻¹ yr⁻¹, respectively). The evolution of erosivity extremes (bi-decadal 98th percentiles, representing an approximation of the hydrological hazard) shows a recrudescence of the hydrological forcing, which started from the change-point year 1490 (after some peaks just prior to 1200), continued throughout the LIA until the end of the 18th century, and then manifested a statistically significant upward trend (Mann-Kendall test p < 0.01) (Figure 7B).

We have analyzed the erosivity density (ED, MJ hm⁻² h⁻¹), which measures the erosivity per rainfall unit and is a better indicator of the erosive hazard than rainfall erosivity (Glur et al., 2013; Diodato et al., 2019a). We addressed this issue by examining the estimated ED time-series since 1891, as it is only from this year that high-quality, well-documented, and homogeneous annual rainfall records are available for the ARB. Figure 8A shows that the erosive hazard has increased significantly from the end of the 19th century until 2016 (Mann-Kendall test p < 0.05). This means that a gradual increase in erosivity density due to extreme rainfall events has increased pressure on the landscape. Although the basin has received a fairly constant amount of precipitation from year to year (Mann-Kendall test p > 0.05), the aggressiveness of rainfall has become even more variable, with extreme values persisting and evolving more unexpectedly. This is in line with the results of Bezak et al. (2020), who have observed an increasing trend in the frequency of peak erosivity years in recent decades (1961–2018) in Italy, particularly in Tuscany. Though not statistically significant, the decadal trend was also positive (Figure 8B).

Evidence for an attribution to climate change of extreme erosive events associated with floods and torrential rains remains limited, partly due to the high spatial variability and unpredictability of these events in the long term. According to Cislaghi et al. (2005) and Pavan et al. (2008), the frequency of rainy days has decreased over Italy over the 19th and especially the 20th centuries, but short-duration episodes (i.e., from one to three hours) have instead enhanced the torrential character of seasonal rains. Colarieti Tosti (2014) also pointed out that in the coming decades the polar vortex is likely to undergo a phase of expansion toward southern latitudes, leading to an intensification of the hydrological cycle in the Mediterranean region. It is also expected from global model projections that the frequency of extreme rainfall events over hazard-exposed landscapes in much of Europe would increase toward the end of the 21st century (Myhre et al., 2019; Rädler et al., 2019). These results seem to reflect the observed increased hazard associated with extreme rainfall in the Mediterranean region, despite a decrease in total precipitation amounts, as claimed by Alpert et al. (2002), Paxian et al. (2015), Drobinski et al. (2018) and Bevacqua et al. (2019), but questioned by Mariani et al., (2014). De Luca and Galasso (2018) also revealed a substantial stability in the temporal and spatial behavior of intense rain events during the 20th century. These contrasting results underline that different analyses can change the perception of the type of hazard associated with hydrological processes when dissimilar metrics are used. Reconstructed rainfall erosivity in two Mediterranean fluvial basins, the Calore Basin in southern Italy (Diodato et al., 2008) and the Po Basin in northern Italy (Diodato et al., 2020c), showed similar increasing trends since the end of the LIA. If this rainfall regime continues, it could result in an increased erosive hazard affecting Mediterranean lands because erosive events occur in a more erratic way. This implies that a future increase in extreme rainfall could have serious consequences in terms of soil erosion, sudden flooding and various types of ecological disruption in the Mediterranean region (Borrelli et al., 2017). An increase of the geo-hydrological hazard due to more intense rainfall activity over the past five decades (1970–2019) is indeed recognized for hilly inland areas of central Italy (Paliaga et al., 2020). However, a sustained long-term trend, dating back to 1760, has not been observed in northern Italy, despite the marked inter-annual and inter-decadal variability of storms (Pieri et al., 2016). This underlines the need for high-resolution climate model projections to reliably predict erosive precipitation. It remains however difficult to determine how Mediterranean cyclones produce erosive trends due to extreme rainfall events (Lionello et al., 2016). Deepening in warm waters, Mediterranean cyclones form mainly around a few centers, with a dominating area in the Gulf of Genoa (∼240 km², centred on 44°10′N, 08°55 E), where a slow-moving low-pressure field (or Vb-weather pattern) can transport large amounts of rain (Hofstätter et al., 2016). While the incidence of cyclones has decreased in the Mediterranean, the increase in erosivity associated with short-term and extreme rainfall could be due to a sharp increase in convective rainfall associated with extreme events during the recent warming. Berg et al. (2013) have shown that convective rainfall is more responsive to temperature increases than frontal passages, and increasingly drives extreme erosive events. Aggressive precipitation is likely to become more frequent as global warming increases (AghaKouchak et al., 2018), albeit with different patterns of change in small areas without the emergence of distinct homogeneous spatial and temporal trends (Libertino et al., 2019).

The wavelet power spectrum (Figure 7C) reveals significant high-frequency periodicities in the erosivity time-series that tend to mirror the well-known ∼11- and ∼22-year sunspot cycles (with ∼22-years magnetic solar cycle consisting of two ∼11-year cycles with opposite polarities; Schwabe, 1843; Hale et al., 1919), and the quasi-secular (∼80–100 years) Gleissberg (1958) cycle. The Sun can considerably affect the hydrological cycle through multi-scale feedbacks possibly originating complex geographical distributions of solar-related signals in hydroclimatic factors. Climate sensitivity to solar variations (magnetic activity of sunspots and the Sun’s orbital dynamics), which follows a frequency-dependent solar energy transfer function (Scafetta and West, 2006), can cause regional changes in the moisture-holding capacity of the atmosphere and alter the hydrological cycle and precipitation patterns (North, 2004). In river basins in particular, the integrated nature and the inertia of river and sediment discharges may reveal solar-induced pulses that amplify precipitation signals, while the links between solar activity and precipitation may be comparably weak (e.g., Milly and Wetherald, 2002; Tomasino et al., 2004).

Our reconstructed storms (translated into rainfall erosivity data) can be considered as suitable proxies of fluvial discharges. Historical sources provide indeed a deal of evidence for the identification of hydrological extremes in the ARB, which includes the discharge of water and sediment from the river. In his Codex Atlanticus Leonardo da Vinci wrote for instance that “the tempest sent a sudden flood of water to submerge al the low part of this city; added to which there came a sudden rain, or rather a ruinous torrent and flood of water, sand, mood and stones, entangled with roots, and stems and fragments of various trees” (Pfister et al., 2009, p. 74). It dates to 1503 the information on “the irregularity in the flow of the Arno waters and their tendency to continuously change the river bed by transporting new sediments” (Curti, 1979, p. 40). Hints for a possible influence of solar cycles on the flood discharges of the Arno river are provided by Zanchettin et al. (2008), who found solar-type periodicities related to the magnetic activity of sunspots and suggested that Sun could be one of the drivers of hydrological processes in northern Italy. In northwest Italy (1701–2019 CE), Diodato et al. (2021) also found statistical links between rain-hazard metrics (erosivity density and return periods of erosivity peak values) and the ∼22-year solar cycle together with Atlantic teleconnection patterns. These studies indicate that regional patterns of precipitation or temperature changes can be modified by large-scale climate structures (teleconnections) induced by variations in solar energy absorption in the atmosphere and ocean surface layers (Perry, 2007). In Europe, the frequency of regional precipitation oscillation peaks can be explained by persistence in atmospheric and circulation patterns over the North Atlantic like sea-surface temperature (SST), as reflected by the Atlantic Multidecadal Oscillation (AMO) index (Willems, 2013). A possible climate mechanism governing the frequency of storm erosivity in the ARB is the Atlantic Multidecadal Variability (AMV), i.e., the change in SST in the North Atlantic. Based on terrestrial proxy records from the circum-North Atlantic region, the AMV reconstruction by Wang et al. (2017) exhibits the pronounced variability on multidecadal time-scales commonly referred to as the AMO (Kerr, 2000; Enfield et al., 2001), i.e. the internally generated component of AMV. We observe that during a warm-dominating phase of the AMV, the number of storms is significantly lower whereas during a cold-dominating phase of the AMV the number of storms is significantly higher (Figure 7D). In this sense, the strength of solar activity apparently modulates the connection between the AMV and the storms in the ARB.

While Sutton and Dong (2012) have identified links between the increase in mean seasonal precipitation in the Mediterranean region and the negative phases of the AMO, it is known that the AMO tends to be related to the North Atlantic Oscillation (NAO) index (Hurrell, 1995; Jones et al., 1997; Cook et al., 2019), i.e. the difference in normalised sea-level atmospheric pressures between stations near Azores (roughly 38°N) and near Iceland (roughly 65°N). Walter and Graf (2002), in particular, have shown that the AMO is negatively correlated with the NAO during its cool (negative) phases while the relationship is weak with the warm (positive) AMO phases. NAO-associated flow moisture convergence anomalies explain part of the precipitation patterns of the Mediterranean region (Seager et al., 2020). Casanueva et al. (2014) examined several Northern Hemisphere teleconnection patterns and found an increase in extreme precipitation events in the central and eastern Mediterranean regions with positive summer values of the NAO index for the period 1950–2010. In a multi-centennial perspective (1680–2019 CE), Diodato et al. (2020a) showed that substantially neutral (winter) NAO states (without a clear dominance of positive or negative phases) correspond to a decline of erosivity extremes throughout the Mediterranean. In the present basin-wide study, we observed (Figure 7E) a predominantly weak anti-correlation between the reconstructed erosivity time-series and the proxy-based multi-annual NAO reconstruction of Hernández et al. (2020). These studies highlight the need for a more detailed understanding of the linkages between large-scale atmospheric and oceanic oscillations and different erosivity outputs (e.g., average or extreme events), which depend on the geographic context (e.g., large regions or small river basins) and the temporal frame (e.g., seasonal or yearly) in which the hydrological response is assessed.