# Millennium-scale changes in the Atlantic Multidecadal Oscillation influenced groundwater recharge rates in Italy

Nazzareno Diodato, Gianni Bellocchi

► To cite this version:

Nazzareno Diodato, Gianni Bellocchi. Millennium-scale changes in the Atlantic Multidecadal Oscillation influenced groundwater recharge rates in Italy. Communications Earth & Environment, 2024, 5 (1), pp.56. 10.1038/s43247-024-01229-6. hal-04955154

HAL Id: hal-04955154

https://hal.inrae.fr/hal-04955154v1

Submitted on 18 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Millennium-scale changes in the Atlantic Multidecadal Oscillation influenced groundwater recharge rates in Italy

Nazzareno Diodato<fnref n="1" /> & Gianni Bellocchi id 1,2✉

Groundwater, essential for water availability, sanitation and achieving Sustainable Development Goals, is shaped by climate dynamics and complex hydrogeological structures. Here, we provide a time-series of groundwater recharge from 801 CE to the present day in the Tiber River Basin, Italy, using historical records and hydrological modelling. Groundwater drought occurred in 36% of Medieval Climatic Anomaly (801-1249) years, 12% of Little Ice Age (1250-1849) years and 26% of Modern Warming Period (1850-2020) years. Importantly, a predominant warm phase of the Atlantic Multidecadal Oscillation, aligned with solar maxima, coincided with prolonged dry spells during both the medieval and modern periods, inducing a reduction in recharge rates due to hydrological memory effects. This study enhances understanding of climate-water interactions, offering a comprehensive view of groundwater dynamics in central Mediterranean and highlighting the importance of the past for sustainable future strategies. Leveraging this understanding can address water scarcity and enhance basin resilience.

<footnote n="1">¹Met European Research Observatory - International Affiliates Program of the University Corporation for Atmospheric Research, Benevento, Italy.</footnote>

<footnote n="2">²Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France. ✉email: gianni.bellocchi@inrae.fr</footnote>
Understanding the importance of groundwater is critical to effectively advancing the Sustainable Development Goals. Groundwater plays a crucial role in sustaining aquatic ecosystems and has gained importance for human development¹,². However, recent decades have witnessed the dual impact of recurring heat-induced droughts in mid-latitude ecosystems³ and unpredictable precipitation trends⁴, resulting in substantial anomalies in current groundwater recharge patterns⁵. While the literature suggests a potential link between these recharge anomalies and the depletion of aquifer storage in various regions worldwide⁶, factors such as decreased infiltration rates and increased evaporation may contribute to the complex interactions between groundwater recharge patterns and aquifer storage dynamics⁷. Thus, the imperative lies in enhancing our understanding of water supply indicators within river basin ecosystems. These indicators not only hold importance for the scientific community by providing guidelines for landscape conservation planning⁸ but also offer crucial insights to governments for evidence-based decisions supported by water supply data⁹⁻¹¹. However, the study of groundwater recharge (GWR) has received less attention compared to surface water hydrology¹², with its exploration predominantly limited to short time scales¹³. This is primarily due to the challenge of obtaining reliable estimates of recharge rates, which remains a persistent issue in various research endeavours¹⁴. Some studies suggest that altered precipitation regimes resulting from changes in land use/land cover can exert substantial impacts on GWR¹⁵, thereby influencing water resource availability and the uninterrupted provision of ecosystem services¹⁶. To illustrate the influence of soil layer positioning on rainwater absorption, the Italian Renaissance polymath Leonardo da Vinci (1452–1510) recognised gravity as the driving force behind water flow on mountain slopes, as captured in the *Codex Parix G* (Pfister et al., p. 59)¹⁷:

> And if ..., the rain, which penetrates from the summit of the mountain to the base, which is the boundary of the sea, descends and softens the slope opposite to the said mountain and constantly draws the water ...

This insight into mountains as Earth's water reservoirs suggests da Vinci's speculation about their potential to supply a considerable portion of essential water to both natural ecosystems and human populations.

The importance of precipitation as a fundamental source for replenishing global water catchments is well-established and reinforced by numerous studies¹⁸⁻²⁰, including research conducted in Italy²¹,²². Notably, the Tiber River, with its indelible mark on the landscape, stands intricately linked to Italy's geography and history²³. Serving as a maritime gateway to the Mediterranean Sea, the river historically acted as a conduit for goods to and from Rome, the hub of the Roman Empire. Similar to other regions around the world, Italy, particularly the Tiber River Basin (TRB), has experienced the consequences of climate change²⁴. This is evident in shifts in precipitation timing and distribution²⁵,²⁶, leading to large interannual water balance fluctuations²⁷,²⁸. This phenomenon is accentuated in mountainous river basins worldwide, where irregular recharge modulation is predominantly driven by precipitation²⁹ (Fig. 1).

Remarkably, while different studies recognise the Mediterranean region's drought as a potential precursor to drought occurrences in other European areas, the decadal variability in the central Mediterranean have received comparatively less attention³⁰,³¹. Although simulating a river basin ecosystem using an integrated and distributed approach, accurately scaling precipitation inputs remains a source of uncertainty, potentially yielding misleading groundwater response outcomes³². Likewise, projections originating from distinct general circulation patterns exhibit notable divergence³³, introducing

uncertainties of greater magnitude than those stemming from multi-parameter ensemble approaches³⁴.

Particularly in regions like South Africa, the Amazon and the Mediterranean basin, discernible drying trends are evident in the CMIP6 projections of Global Climate Models. However, the extent of this drying remains contingent on climate model outcomes³⁵. On a smaller scale, historical hydroclimatic variations, particularly the intricate interplay between climate and groundwater dynamic, remain insufficiently understood³⁶,³⁷. The complexities inherent in reconstructing long-term GWR data further contribute to this limited understanding³⁸. As such, the historical reconstruction and prediction of GWR data continue to present challenges across a range of temporal scales, encompassing existing Earth systems³⁹, environmental engineering and hydrogeology fields⁴⁰,⁴¹. Consequently, the resulting variations in groundwater feedbacks to the climate system remain inadequately understood⁴², impeding the establishment of meaningful interdisciplinary connections.

For better future planning, a more nuanced comprehension of how groundwater recharge has responded to historical climate shifts is essential, transcending solely human-induced impacts¹³. Within the hydrological cycle, climate change's impact extends to rainwater infiltration, deeper percolation, and consequently, groundwater recharge, operating on timescales from decades to millennia (Fig. 2e). In the pursuit of advancing groundwater recharge studies, this study introduces a distinct approach and methodology that bridges the gap between insights within water research and limited available data for building a historical model. To achieve this, we present the Historical Hydrological Recharge Model (HHydroReM), a parsimonious model designed to reconstruct the annual groundwater recharge rate (GWR) specifically for the TRB region. This model not only addresses data scarcity challenges but also contributes to a more comprehensive understanding of groundwater dynamics.

The Tiber River Basin holds considerable hydrological importance within the Italian peninsula, making it an ideal case study for gaining insights into water resources in the region (Fig. 2a, b). The Italian peninsula's diverse geography and climate, stretching from the northern Alps to the southern Mediterranean coastline, add complexity to its water dynamics. As a major watercourse in central Italy, the TRB encompasses regions like Tuscany, Umbria and Lazio, playing a crucial role in supplying water for agriculture, industry and human settlements. The study of TRB offers key insights into hydrological processes, groundwater dynamics and the intricate interplay between surface water and groundwater. This understanding is essential for ensuring sustainable water availability and stability throughout the Italian peninsula. The TRB serves as a representative model for water resource management challenges prevalent in the region. Its hydrological conditions range from the relatively wet northern regions to the drier southern areas, rendering it an invaluable case study for assessing climate variability's impacts on groundwater recharge rates, water storage and overall availability. Moreover, the Italian peninsula's central location within the Mediterranean Sea renders it representative of the broader central Mediterranean region. The peninsula shares environmental and hydrological characteristics with the Mediterranean basin. The cyclical nature of rainfall, leading to water surplus during the wet season (Fig. 2d), has shaped the millennia-old river valley of the Tiber (Fig. 3a, b), reflecting the distinct features of the Mediterranean region. Specifically, the hydrological dynamics of central Italy's GWR are substantially influenced by abundant precipitation patterns. The TRB receives an average annual precipitation of about 1000 mm (1981–2010). However, certain parts within the southeastern and northern areas can encounter even more substantial levels of precipitation, occasionally surpassing 1200 mm yr⁻¹ (Fig. 2c).
<figure><img src="image_3.png" /></figure>

**Fig. 1** Overview of global groundwater resources. This map offers insights into the abundance, quality and vulnerability of Earth's groundwater resources, serving as a tool to convey groundwater-related concerns to hydrology experts, decision-makers and the general public. The central Mediterranean region receives substantial recharge in local aquifers with complex hydrogeological structures, demanding a comprehensive understanding, knowledge and careful planning and management due to limited available information regarding this hidden resource in many areas. Image adapted from the World-wide Hydrogeological Mapping and Assessment Programme of the German Federal Institute for Geosciences and Natural Resources and the United Nations Educational, Scientific and Cultural Organization; image source: https://www.whymap.org/whymap/EN/Home/whymap_node.html ("This information may be freely used and copied for educational and other non-commercial purposes, provided that any reproduction of data be accompanied by an acknowledgement of the source", https://www.whymap.org/whymap/EN/Disclaimer/disclaimer_node_en.html;jsessionid=515B048B1383C0E75A79F48AEAB62849.internet992).

The monthly regime exhibits a dry phase during the summer, resulting in a partial water deficit (Fig. 2d, represented by beige and brown colours). This is succeeded by a period of positive change in the water deficit, marked by soil water replenishment (Fig. 2d, depicted in blue light). The subsequent extended winter phase experiences a water surplus (shown in blue), which facilitates infiltrations, consequently recharging the groundwater of the TRB. This cyclic water balance pattern is a consequence of the Mediterranean climate, alternating between quiet summer conditions occasionally punctuated by short stormy episodes, and longer rainy periods in winter and spring, characterised by water surplus and runoff.

The HHydroREM model adopts a deliberately sober approach, relying on simplified inputs such as intensity indicators, duration and inter-storm drought intervals⁴³. This design choice aims to overcome the constraints posed by more complex, physically-based models. The model's estimations were substantiated through extensive validation, leveraging hydrological measurements of groundwater capacity from two TRB sites: Peschiera and Boschetti (Fig. 2c), spanning the period 1941–1977. Following calibration and validation, the model was employed to assess the TRB's annual recharge across distinct climatic periods: the Medieval Climatic Anomaly (MCA, 801–1249), the Little Ice Age (LIA, 1250–1849) and the Modern Warming Period (MWP, 1850–2020). This endeavour has yielded the longest time-series of annual GWR data (801–2020 CE), covering a geographical scope suited for hydrological studies in the central Mediterranean and beyond. The TRB's role as a case study holds promise for shedding light on hydrological responses to climate fluctuations and the intricate water resource management challenges faced by the broader central Mediterranean region.

## Results and discussion

**Model calibration and validation.** We calibrated our HHydroREM model using recharge data from the modified Thornthwaite groundwater recharge water balance for the period 1928 to 1974. This 47-year calibration dataset (y) revealed a significant linear relationship (y = a + b·x) with the predicted (x) dataset (F-test, *p* ~ 0.00), with a high R² value (0.74). The parameters of Eq. (2) were determined as follows: A = 9.137 (a parameter of scale), converting the term within brackets to mm yr⁻¹ and B = 202 mm yr⁻¹, representing the millennium baseline flow (Fig. 2f). Additionally, η = 2.00 is an empirical coefficient, while α = 4.00, β = 21.00 and γ = 0.50 are parameters modulating storm precipitation intensity during recharge months. The intercept (a = −0.719 ± 41 standard error) and slope (b = 1.000 ± 0.087 standard error) yielded a single data-point (1929) outside the 95% prediction bounds for new observations (Fig. 4a, light pink band). The mean absolute error of 68 mm yr⁻¹ was smaller than the standard error of the estimates (78 mm yr⁻¹).

We assessed the normality of the GWR(H)–model data distribution, represented by Eq. (2), using standardised skewness and kurtosis (Table 1). Deviations beyond the −2 to +2 range indicate significant departures from normality. Both standardised skewness and kurtosis values for the calibration dataset remained within the acceptable range, reflecting no issues with the distribution of the calibration data. While the standardised kurtosis values for the GWR(H)–model aligned with expectations, indicating acceptable tail behaviour, the standardised skewness value exceeds the normal range, indicating asymmetry in data distribution. While we should be mindful of this asymmetry when interpreting results from the GWR(H)–model, the Kolmogorov–Smirnov (K-S) test indicates that the two data
<figure><img src="image_4.png" /></figure>

**Fig. 2** Environmental, climatic and hydrogeological features of the study area. **a** Mediterranean region framed within the global context. **b** Detailed view of the Mediterranean area with emphasis on the central Italian peninsula, highlighted by the red box. Image source for **a** and **b**: ClipSafari (https://www.clipsafari.com), available under the Creative Commons Attribution 0 license, https://www.clipsafari.com/clips/0188901-a-planetary-perspective-a-view-of-the-earth-from-space. **c** Distribution of mean annual precipitation (1971-2000) across central Italy, featuring the delineation of the Tiber River Basin (highlighted by the red square). Location of monitoring stations, Boschetti well groundwater and Peschiera discharge, are also marked. Image generated using Climate Explorer, https://climexp.knmi.nl, with data from GPCC-Global Precipitation Climatology Centre Full Data Monthly Product Version 2020 for 0.25° resolution (https://opendata.dwd.de/climate_environment/GPCC/html/fulldata-monthly_v2020_doi_download.html). **d** Representation of the Thornthwaite modified hydrological water balance over the Tiber River Basin during 1981-2010. The calculation is based on GPCC-Global Precipitation Climatology Centre (https://climatedataguide.ucar.edu/climate-data/gpcc-global-precipitation-climatology-centre) and CRU-Climate Research Unit of the University of East Anglia temperature (https:// rudata.uea.ac.uk/cru/data/hrg/cru_ts_4.04) data. Different shades illustrate water soil deficit (brown), negative change in soil moisture (beige), positive change in soil moisture (blue light), water surplus (blue) and precipitation (blue). Adapted from the WebWIMP Programme - WebWIMP 1.01, implemented at the University of Delaware, USA, in 2003 (source: http://climate.geog.udel.edu/~wimp/index.html, accesses on December 2022). **e** Depiction of the dynamic phases of a typical hydrological cycle (with component labels in orange) and the long-term pattern of groundwater flow over decades to millennia, showcasing aquifers as a 'natural buffer' against drought impacts. The graph was created by the authors, using composite images primarily sourced from a NOAA (National Oceanic and Atmospheric Administration) image (https://www.noaa.gov/education/resource-collections/freshwater/water-cycle - image credit: Dennis Cain/NW) and incorporating elements from an image published by IAH (International Association of Hydrogeologist, Strategic Overview Series: Climate-Change Adaptation & Groundwater) in a freely available report (https://iah.org/wp-content/uploads/2019/07/IAH_Climate-ChangeAdaptationGdwtr.pdf). Information within the NOAA repository is in the public domain (https://repository.library.noaa.gov/Content%20and%20Copyright).
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Historical scene of Tiber valley and its river. **a** View of the mid-19th century Tiber valley from Monte Mario (Rome, Italy) with the Milvio bridge, as portrayed in a painting by the Russian painter Matveev Fëdor Michajlovič (1758-1826). Image source: Koller Auctions, lot 3086/A190, https://www.kollerauktionen.ch/it/object-search.htm?term=FEODOR%20MICHAILIOWITSCH%20MATVEEV&origin=. Permission for the image's use was kindly granted by Laura Järmann of Koller Auctions, Zurich, Switzerland. **b** As **a**, with a focus on the Tiber River course crossing the urban sector of Rome City, 1690s, in a painting of the Dutch painter and draughtsman Gaspar van Wittel (1652/1653-1736). Image source: https://it.wikipedia.org/wiki/Tvere#/media/File:The_Castel_Sant'Angelo_from_the_South.jpg (available under the Creative Commons CC0 License; https://creativecommons.org/publicdomain/zero/1.0/).</figcaption></figure>

datasets (calibration and modelled) likely share the same distribution, supported by a small maximum distance between distributions (Dn = 0.13), a two-sided large sample K-S statistic of 0.62 and a high p-value of 0.84. Figure 4b shows a Gaussian-like distribution of the model residuals, devoid of skewness. Further reinforcing this, the Durbin-Watson (DW) statistic (DW = 2.08) with a p-value exceeding 0.05 (p = 0.60) indicates no serial autocorrelation in residuals. The validation dataset (1999–2020; Fig. 4c) demonstrated a significant linear relationship with the predicted data (F-test, *p* ~ 0.00; *r* = 0.71), and no serial autocorrelation in residuals (DW = 1.53; *p* = 0.13).

**Extended model validation.** The validation was extended to incorporate measurements of groundwater capacity, comparing model estimates with the hydrogeological response of the TRB. We used standardised mean values of the water-table level measurements performed at Boschetti well⁴⁴ and Peschiera spring⁴⁵ in 1941–1974 for this central test of the model’s effectiveness. Indeed, as in the absence of direct observations, recharge data are most often derived from simplified hydrogeological budgets (as in this study) without a priori check of the model accuracy⁴⁰.

To make modelled discharge values statistically comparable to water-table level measurements, we transformed the two variables into two time-series of relative standardised anomalies. Their joint evolution (Fig. 4d) demonstrated a significant linear relationship (F-test, *p* < 0.05, *r* = 0.72) between modelled recharge values and the hydrogeological response of the TRB. The two times-series in Fig. 4d mostly overlap and the bars of the model residuals (Fig. 4e) reflect a normal distribution (blue line). There is not a statistically significant difference between the two distributions (maximum distance between the two distributions Dn = 0.16, two-sided large sample K-S statistic = 0.70, *p* = 0.72). Statistical comparisons and visualisations supported the operational use of the HHydroReM.

**Historical evolution of groundwater recharge.** Throughout history, spring water availability from GWR has been a consistent natural phenomenon for the Tiber basin, constituting an enduring natural occurrence. These preliminary findings suggest that GWR never dropped below a minimum level, set at around 150 mm per year, sustaining springs feeding into the Tiber, which in turn also flow into the river itself, thereby regulating the Tiber continuously. The Italian engineer Elia Lombardini (1794-1878)

gives proof of this in his historical work, Guida allo studio del-l'idrologia fluviale e dell'idraulica pratica (1870), attributing "the uninterrupted Tiber flow to numerous large cavities which are believed to exist in the bowels of the mountains constituting its basin; in which a large part of the rainwater and the waters coming from the liquefaction of the snow would be collected in a large reservoir to flow into the river in the form of springs" (Betocchi, pp. 31-32)⁴⁶. More recently, Vaquero Piñeiro⁴⁷ also discussed the contributing factors to flooding in Rome, attributing a notable role to the inflow of groundwater during flow events. Groundwater drought (GWD) severity is commonly described using a threshold level concept via a deficit index⁴⁸. In this approach, a GWD event begins when the reading falls below a set threshold (e.g., 25th percentile) and ends when it rises again⁴⁹. This percentile is useful to assess the GWD occurring in the three climatic periods: the Medieval Climate Anomaly (MCA, 801-1249 CE), the Little Ice Age (LIA, 1250-1849 CE) and the Modern Warming Period (MWP, 1850-2020 CE). Unfortunately, the scarcity of historical evidence in medieval times frequently makes determining how and why climate and weather affected water supply challenging. To summarise the above three phases, we obtain percentages of years with values below the 25th percentile: 36, 12 and 25%, respectively for the MCA, LIA and MWP. From these percentages, it can be inferred that the GWD that occurred during the LIA period had a lesser impact on the underlying groundwater, compared to the MWP, and secondarily to the MWP. In fact, it is during the medieval period that we must refer to in order to find the most severe and recurrent GWD events that have affected the TRB (Fig. 5d, yellow bands). We recall, for instance, the one that spanned, on several occasions, several decades between the years 880-980 CE and 1009-1160 CE, when all of Europe was involved in a multidecadal GW megadrought⁵⁰,⁵¹. In Italy it is reported towards the end of this period, i.e. in the year 1131, with great heat, drying up of all rivers and famine⁵². The maxima of the medieval period, stands, exceptionally, at 814 and 849 mm, and dates back to the year 860 and 1240 CE, respectively. Figure 6a shows the spatial distribution of this drought over the central Mediterranean area using the Palmer Drought Severity Index (PDSI) from data by Cook et al.⁵³, confirming hydrological drought (indicated by orange to red colour bands) over central and northern Italy, which were more pronounced during the period 1040-1100 CE. Interestingly, such drier climatic conditions are visible also in the isotopic data from Pergusa (Sicily) at around 1050-1150 CE⁵⁴. If we observe the corresponding evolution of the Atlantic Multidecadal Oscillation
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Groundwater recharge model calibration and validation. **a** Scatterplot comparing modelled groundwater recharge and Thornthwaite's groundwater (GW) recharge's water balance using Eq. (2), during the calibration period of 1928–1974, 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). **b** Related histogram of residuals and its normal fit (blue line). **c** Coevolution between the modelled GWR(H) (represented by the orange line in mm) and the groundwater level (indicated by the blue line in m) of the Castelporziano piezometer E13.2 during the period 1999–2020 (data sourced from Banzato et al.⁹⁷); **d** Coevolution between the modelled GWR(H) (orange line) and hydrogeological response (actual discharge) of TRB (blue line) standardised anomalies at the validation results (data sourced from *Annals Project* of ex SIMN, http://www.bio.isprambiente.it/annalipdf) for 37 years (1941–1974). **e** Relative histogram of residuals, with over-imposed theoretical Gaussian shapes (blue line).</figcaption></figure>

<table><caption>**Table 1** Summary statistics for calibration and modelled data.</caption><thead><tr><td>Statistics</td><td>Calibration dataset</td><td>Model estimates, Eq. (2)</td></tr></thead><tbody><tr><td>Count</td><td>47</td><td></td></tr><tr><td>Mean (m)</td><td>451</td><td>452</td></tr><tr><td>Standard deviation (m)</td><td>152</td><td>131</td></tr><tr><td>Coefficient of variation (%)</td><td>0.34</td><td>0.29</td></tr><tr><td>Minimum (m)</td><td>105 (1945)</td><td>228 (1950)</td></tr><tr><td>Maximum (m)</td><td>875 (1960)</td><td>840 (1960)</td></tr><tr><td>Range (m)</td><td>770</td><td>612</td></tr><tr><td>Standardised skewness</td><td>1.75</td><td>2.79</td></tr><tr><td>Standardised kurtosis</td><td>1.20</td><td>1.72</td></tr></tbody></table>

(AMO) throughout the period of the MCA, we can deduce that this index was always above zero and especially in the phases of greatest drought it remained around +0.3, between 950 and 1100 CE (Fig. 5d, orange line).

The AMO refers to a warming/cooling (AMO > 0/AMO < 0) pattern of multidecadal North Atlantic sea-surface temperature (SST) anomalies with a periodicity of ~60-yr periodicity⁵⁵. SST anomalies, as represented by the AMO, also induce atmospheric pressure gradients, redistributing air masses through the North Atlantic Oscillation (NAO) between the subtropical and subpolar latitudes of the North Atlantic, and modulating the strength and latitudinal position of westerly flows⁵⁶,⁵⁷. In particular, large-scale synoptic configurations conducive to precipitation in the central Mediterranean are generally associated with negative phases of the NAO (NAO-), an oscillating pattern of mean sea-level pressure between the subpolar and subtropical latitudes of the North Atlantic⁵⁸. In the late MCA another major water shortage is reported between the years 1246 and 1269, which is also mentioned in climate reconstructions from a tree-ring network across the Mediterranean indicate drought in Italy⁵⁹. In this period, in fact, approximately 58% of the years had GWR values below the 25th percentile across the TRB.

The MCA also coincided with the Great Solar Maximum between 1100 and 1250 CE and its temperature optimum occurred between 1000 and 1249 CE, when a strong temperature

maximum occurred over most of Europe⁶⁰. As observed in Fig. 5b, c, during the MCA and MWP, GWD events are generally associated with warm climatic phases. In contrast, the LIA phase (1250–1849 CE) shows dry periods that are typically linked to colder conditions, with GWD percentages lower (12%) than those during the MCA and MWP. Figure 5b (blue line) also displays the variability of higher recharge values, which generally correspond to humid periods. Thus, maximum GWR values are expected during the LIA, with the absolute maximum of 926 mm, in the year 1590, preceded and followed by same high values as in the year 1566, with the 595 mm, year 1589 with 402 mm, years 1598 and 1599, with 597 and 496 mm, respectively. These wet periods, characterised by peaks in GWR, align with the temporary advances of Lake Trasimeno (43° 08′ N, 12° 06′ E, near Perugia), as recorded in the year 1602⁶¹ (Fig. 5a). These advances were likely triggered by continued storm events⁴³.

While detailed information on vegetation, climate and anthropogenic impacts in Lake Trasimeno underscores the unique contributions and challenges associated with obtaining relevant data for a detailed GWR analysis, it is essential to note that the Trasimeno lake record, being historical, may have limitations⁶². Until the mid-17th century, important GWD are chased in a sporadic way, with the exception of the longest occurring between 1397–1409 and 1481–1492 CE. However, the advent of the Maunder Minimum (MM, 1645–1745 CE)⁶³ marked a turning point. Between 1663 and 1673, a new dry phase emerged. Subsequently, a period of severe and prolonged GWD resumed after the end of the MM (1715) until 1745, with ~52% of the years below the 25th percentile. This trend can also be observed in the respective hydrological drought in Fig. 6b, which was likely among its causes. There are records, testimonies and news from the year 1746 about wells and streams running out of water, trees and grasses drying up⁵². This period coincides with a phase of low water levels in Lake Trasimeno, which began in 1630 and lasted until about the middle of the 18th century⁶⁴. There are also records of the drying up of Lake Fucino (41° 59′ N, 13° 32′ E) in western Abruzzo (central Italy) around 1752⁶⁵ (De Sanctis, p. 6), with continuous invocations for rain to the Madonna di Capodacqua⁶⁶, whose sanctuary stands near the source of Velino (a sub-tributary of the Tiber), in the
<figure><img src="image_7.png" /></figure>

**Fig. 5** Environmental hydrology evolution over the period 801-2020 CE in the Tiber River Basin. **a** (dashed red line) Fucino lake level anomaly from 501 to 1400 CE⁹⁸ and (continue red line) Trasimeno lake anomaly level from 1401 to 2000 CE⁶⁸. **b** (dark blue line) Annual reconstructed GWR (801-2020 CE) with key climatic periods and drought events (vertical yellow bands). **c** Coloured bands representing surface air temperature anomalies relative to the period 1961-1990 climate baseline for Europe⁹⁹. **d** Coevolution of the Atlantic Multidecadal Oscillation (orange line, source: Mann et al.¹⁰⁰), and the 10th percentile of GWR (light blue line). **e** Wavelet spectrum of the GWR time-series with bounded colours (0.05 significance areas) within the bell-shaped, black contour marking the limit between the reliable region and the region below the contour where edge effects occur, a.k.a. the cone of influence. The violet vertical line passing through the year 1590 divides the time-series, representing the lower extreme of GWR increasing and then declining. Simultaneously, we can observe the reverse trend for the AMO.

municipality of Cittareale (42° 37′, 13° 10′). This long phase is also in harmony with what has been described by Ionita et al.⁶⁷, who found longer but less warm droughts for the LIA than modern ones (Fig. 5c). From the end 1500 s (violet vertical line), the lower end of the GWR (10th percentile) started an important and unstoppable decrease (significant Mann-Kendall trend-test, *p* < 0.05) until the end of the time-series in 2020 (Fig. 5d, light blue line).

This decline in GWR corresponds to the evolution of Trasimeno lake level. From 1400 to the beginning of 18th century, there was a period of stability with minimal oscillations, followed by a gradual descent until the mid-18th century. Afterward the level suffered a drastic decreased to get to present day (Fig. 5a, red line). With the exit from the LIA, the GWD get back to enhance, with the percentage of the years below the 25th percentile equal to 26%, though peaks of a certain importance are expected, but only between the second half of modern era and, precisely, in the year 1937, with 831 mm, year 1941, with 701 mm, year 1960, with 875 mm, and with last peak above 600 mm in the year 1984 with 652 mm. During this period, the

Trasimeno lake's level began a series of swings, and at the end of the 1950s, there was a severe water crisis, during which the maximum depth of the lake dropped to less than 3 m: the changing climate and the presence of the new outlet, which prohibits any storing of water above the outlet threshold, has since caused the average lake level to be too low⁶⁸. This drying is reverberated especially GWR at lower end (10th percentile level), undergoing a notable decline after the 1980s (Fig. 5d, light blue line). This depletion was driven, among other factors, by drought mapped in Fig. 6c (with orange and red colours scattered not only in TRB but also in other parts of Italy). However, AMO appears to influence the lower end of GWR (10th percentile level) throughout the entire time-series (Fig. 5d, orange and light blue lines, respectively). It is essential to note that this influence is not observed on a year-by-year basis, but rather manifests as a coevolution at the decadal and multidecadal levels. Due to various factors, particularly the delayed response of GWR following the end of the LIA, there seems to be a lag, likely associated with decreased flow rates and the input of water from snowmelt. This delay becomes more pronounced post-LIA due to the lowering of the water table, which
<figure><img src="image_8.png" /></figure>

**Fig. 6** Spatial distribution of the median value of the Palmer Drought Severity Index (PDSI) across Italy during the periods of major groundwater recharge depletion. a 1085–1140 CE. b 1716–1769. c 1988–2020 (generated using Climate Explorer, https://climexp.knmi.nl, with data from Cook et al.⁵³, and updated CRU - PDSI 4.05early data). Orange to red colours indicate increasing drought. The PDSI is a widely used drought index that incorporates both precipitation and temperature data into a simplified realistic water balance, considering both water supply (rainfall or snow water equivalent) and demand (temperature, transformed into evapotranspiration units). This index varies between -4 (extreme drought) and +4 (extremely wet), and provides insights into drought severity. The grey rectangle indicates the TRB.

decelerates the water cycle. These declines align with a broader, long-
term decrease in river discharge observed over the past five centuries
in the Mediterranean region⁶⁹. Given the distinctive trend observed in
the lower end (10th percentile level) throughout the historical time-
series and its association with Lake Trasimeno, it is interesting to
investigate the extent to which interannual recharge persists in the
groundwater memory.

To investigate this issue, we estimated the Hurst (H) exponents
of the 10th percentile time-series using the R/S method⁷⁰ and the
residual variance-ratio method⁷¹. The resulting H values were
determined as 0.88 and 1.11, respectively. These values indicate a
predictable time-series because H > 0.65⁷², pointing to a connec-
tion between the time-series and long-term memory. This
connection reflects cyclical influences on the occurrence of
precipitation, driven by large-scale weather systems.

This not only reaffirms the intimate relationship between surface
and subsurface hydrology but also highlights a form of memory in the
propagation of the climatic signal within the groundwater dynamics
across the TRB. Moving on to the analysis, we examined the wavelet
power spectrum (Fig. 5e), which reveals both high-frequency (~22
years) and low-frequency (~300 years) periodicities (starting from the
LIA period onwards). This analysis provides insights into the cyclical
patterns and variations that have influenced the groundwater
dynamics in the central Mediterranean region. Indeed, the ~22-year
sunspot cycle, mirroring the ~22-year magnetic polarity cycle of
sunspot activity, consisting of two consecutive ~11-year sunspot
cycles⁷³, suggests that the solar activity may be among the precursors
of hydrological processes in the TRB, potentially related to the
magnetic activity of sunspots (as similarly identified by Zanchettin
et al.⁷⁴ in northern Italy). Additionally, the ~300-year cycle may
correspond to a quasi-periodic feature of sunspot variability, which
likely influences changes in precipitation regimes in the North Atlantic
region (as inferred by Ojala et al.⁷⁵ from paleo-proxy records
in Finland).

Conclusions

This research highlights the synergy between historical hydro-
logical knowledge and modelling techniques, bridging past, pre-
sent and future hydroclimatic conditions. The HHydroReM
model empowers water managers to mitigate risks, enhance
resilience and ensure sustainable water resource management.

This integrated approach supports current decision-making and
lays a foundation for adaptive strategies to address evolving
hydrological challenges. Notably, the successful evaluation of the
HHydroReM model against historical groundwater recharge data
in the Tiber River Basin holds substantial implications for
effective water management and strategic planning. By robustly
replicating past recharge patterns, the model becomes a valuable
tool for water resource managers, urban planners and policy-
makers, helping to track progress towards the Sustainable
Development Goals (SDGs) and advancing us toward a sustain-
able future. Combining historical hydrological insights with
contemporary computational capabilities, this research bridges
traditional knowledge and modern tools, fostering well-informed
decision-making. This connection offers a comprehensive
understanding of complex hydro-meteorological interactions,
climate oscillations and groundwater recharge, directly con-
tributing to SDGs related to clean water and sanitation (goal 6),
sustainable cities and communities (goal 11), and climate action
isolated (goal 13).

The practical implications of this research radiate across water
resource management stakeholders. Water utilities stand to
optimise infrastructure investments, anticipate water availability
and design drought management plans using a deeper under-
standing of historical recharge fluctuations. Urban planners
can formulate sustainable growth strategies that account
for future groundwater availability changes, aligning with the
inclusive, safe and resilient cities goal (goal 11). Agricultural
sectors can refine irrigation practices and reduce environmental
impact.

Furthermore, the study's emphasis on long-term memory within
groundwater dynamics highlights the need for adaptive management
strategies, answering the call to deal with climate change in goal 13.
Recognising the cyclical influences of sunspot activity and climatic
phases, decision-makers can pro-actively prepare for potential
hydrological pattern shifts affecting water supply. While the applic-
ability of our model to various study areas may be influenced by the
availability of observational data, and future efforts are needed to
explore ways to enhance the model's adaptability to diverse settings
with varying data richness, we highlight the HHydroReM model's
ability to simulate historical droughts and periods of abundance aids
in anticipating similar trends, offering valuable foresights.
## Methods

**Data source and variables.** The Tiber River Basin (TRB), spanning over 17,000 km² in central Italy, was chosen as the study area. The TRB is home of the third longest river in Italy (the longest in central Italy), extending for 405 km through the heart of the country. Altitudinal variations within the TRB range between 400 and 1000 m a.s.l., encompassing over 40% of the basin. The basin's topography features slopes predominantly falling between 20% and 40%, rendering it a well-suited candidate for thorough analysis. A consistent set of annual areal mean precipitation and annual mean temperature data for the TRB during the period 1928–1974 were obtained, enriched by incorporating surface runoff data sourced from the *Annals Project* (http://www.bio.isprambiente.it/annalipdf), originally overseen by the former Italian National Hydrological and Oceanographic Service (SIMN – *Servizio Idrografico and Mareografico Nazionale Italiano*). This dataset was updated and made available by the Italian National Institute for Environmental Protection and Research (ISPRA – *Istituto Superiore per la Protezione e la Ricerca Ambientale*; https://www.isprambiente.gov.it). The characterization of soil-moisture storage capacity was drawn from statistical records⁷⁶.

**Thornthwaite's groundwater recharge's water balance.** Groundwater recharge refers to water added to an aquifer through the unsaturated zone after infiltration and percolation following any rainfall event⁷⁷. GWR rates are the rates at which infiltrating water moves across the water table. At a given site, GWR rates can be determined using different approaches, including hydrogeological water balance, one-dimensional soil water flow, groundwater level fluctuation, groundwater balance, and isotope and solute profile techniques⁷⁸,⁷⁹. Thornthwaite's water balance model was used to estimate GWR groundwater recharge for the TRB. This model, depicted in Fig. 7, evaluates the various climatic and hydrological inputs based on Thornthwaite's framework⁸⁰⁻⁸². Notably, Durga Rao et al.²⁹, Vandewiele et al.⁸³ and Mammoliti et al.⁸⁴ have demonstrated that, despite its poor data requirements, the Thornthwaite-revised approach remains a practical and accurate method for estimating monthly water balance in ungauged catchments. The model inputs include mean areal precipitation (mm), mean areal temperature (°C), direct runoff factor from surface runoff and soil-moisture storage capacity. Location-specific

parameters were considered along with additional adjustable parameters (Fig. 7). Thornthwaite's primary emphasis is on assessing water surplus. A detailed breakdown of the water balance components can be found in McCabe and Markstrom⁸⁵. The water balance tool can be freely downloaded from https://pubs.er.usgs.gov/publication/ofr20071088.

Before proceeding with the calibration of the historical GWR model, a crucial step involves verifying the consistency between the recharge assessed through the Thornthwaite-revised approach and the estimated areal basin-scale recharge of the Tiber River aquifer. To achieve this, we conducted a comparison between the Thornthwaite-derived output and the water equivalent thickness data across TRB from the GRACE-NASA MEaSUREs Program⁸⁶, which provides a reliable indicator of terrestrial water storage estimates. In Fig. 8a, this comparison is depicted, with the Thornthwaite-derived output represented by the light blue histogram and the water equivalent thickness shown as red dots. The data span the period from 2004 to 2018. The resulting correlation coefficient is 0.85, meaning high significance (F-test, *p* < 0.05). This is evident from the scatterplot shown in Fig. 8b. Additionally, normalised data were used to illustrate the QQ-Plot in Fig. 8c, which shows that the Thornthwaite-based output data follow a normal distribution, as indicated by the alignment of dots along the relative bisector. Based on these findings, we have reasonable grounds to assume that the interannual changes in the aquifer are reflected in the variability of the Thornthwaite-derived output. This, in turn, enables the calibration of the historical model GWR(H), with a highly noteworthy level of approximation.

**HHydroREM historical model.** The sequence of processes that summarises the GWR propagation denotes the change of the recharge signal as it moves through the terrestrial phase of the hydrological cycle⁸⁷. During dry spells, drainage and runoff levels tend to be low. However, potential evapotranspiration may intensify, leading to increased actual evapotranspiration and an additional loss of water from soil and open water bodies. Eventually, soil moisture storage depletion reduces recharge into the aquifer system, resulting in a reduction of groundwater levels. Actual groundwater levels are influenced by pre-event circumstances and the rate of decline, which, in turn, is determined by the volume of recharge and outflow, as well as the aquifer's storage properties. As groundwater's response to climatic input often exhibits delay and smoothing, groundwater anomalies frequently manifest the results as extended periods of above or below-normal hydrological and climatological variables. The historical estimation of groundwater recharge (mm yr⁻¹), denoted as GWR(H), was performed using a non-linear multivariate regression model tailored to the TRB. This concept was derived from Oke⁸⁸, who established a relationship between GWR and corresponding monsoon-season storm values using a non-linear regression equation. As a result, even in years with low total precipitation, storm precipitation may contribute to some of the highest recharge, as seen in both arid and wet regions across the world⁸⁹. This model incorporates a range of hydro-meteorological factors, including the monthly storm severity index (MSSI), snow severity index (SSI) and Palmer Drought Severity Index (PDSI) data. The formulation of GWR(H) followed a parsimonious approach, integrating both current-year hydro-meteorological conditions and climatic factors over multiple years. These data serve as proxies that, upon processing, can be correlated with Thornthwaite's data and used for estimations in historical times, such as centennial to millennial epochs, for which monthly precipitation and temperatures are not available as inputs for the Thornthwaite model.

<figure><img src="image_9.png" /><figcaption>**Fig. 7 Schematic representation of the water-balance model.** The model was assessed using the Graphical User Interface (GUI) module of the modified Thornthwaite model (adapted from McCabe and Markstrom⁸⁵).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8 Groundwater recharge dynamics: Thornthwaite-based versus GRACE analysis. a Coevolution of groundwater recharge using the Thornthwaite-revised approach (depicted by the light blue histogram) and water equivalent thickness (represented by red dots). b Scatterplot between normalised water equivalent GRACE and normalised Thornthwaite's groundwater recharge water balance using - Eq. (2). c Distribution of sample versus theoretical quantiles of the QQ-plot of Thornthwaite-based output data.</figcaption></figure>

Based on the seminal function of Diodato and Fiorillo⁹⁰, it includes two independent variables, X and Y, representing respectively a non-linear hydro-meteorological function, f(HMF) for the current year (y = 0), and a climatological factor, ClimF, which includes more years, expressed as follows:

$$GWR(H)_{y=0} = A \cdot [f(HMF)_{y=0} + (\text{climF})] + B \quad (1)$$

where GWR(H)_{y=0} is the groundwater recharge (mm).

Expanding the two main terms, we obtain:

$$GWR(H)_{y=0} = A \cdot \left\{ \begin{array}{l} \left[ \text{MSSI}_{j} \cdot \sum_{j=8}^{j=4} \left( 1 + \gamma \cdot \cos\left(6.28 \cdot \left(\frac{j+a}{b-j}\right)^{\eta}\right) \right) \right]_{y=0} \\ + \eta \cdot \sum_{y=11,j=8}^{y=10} \text{MSSI}_{y,j} + \sum_{y=2}^{y=0} \text{SSI}_{y} + \sum_{y=1}^{y=0} \text{PDSI}_{y} \end{array} \right\} + B \quad (2)$$

The term enclosed within the square brackets denotes the non-linear hydro-meteorological function. Subsequent terms indicate the climatic influences, consisting of three components that drive the recharge process. In particular, the MSSI as outlined by Diodato et al.⁴³, incorporates both precipitation intensity and duration. It accounts for the total precipitation from August to April in the current year (y = 0) and is further weighted by a soil moisture index (indicated by the term in round brackets). The succeeding MSSI term operates over the previous 11 years leading up to the current year (y = -11 to y = -1), specifically for the months of August, September and October. This temporal arrangement identifies a memory effect centred on the current year (y = 0), facilitated by an 11-year moving average window function. Our approach draws inspiration from Esit et al.⁹¹, who highlighted the role of multi-seasonal to multi-year memory in deep soil and groundwater moisture for enhancing hydro-climate predictability. Notably, we assign maximum weight to the year y = 0 to robustly account for interannual variability (reflected in the first term of the numerator). In this way, storm-induced precipitation can contribute substantially to recharge, even in years with lower overall precipitation. This phenomenon has been observed in both dry and wet regions⁸⁹,⁹². The second term in the numerator captures the positive shift in soil water deficit as the transition occurs from dry to wet season in September-October. To achieve this, the quantity of rainfall received by the basin during these autumn months becomes of paramount importance. During this time of year, infiltration capacity reaches its maximum potential due to the drying of the soil during the preceding summer months. This drying process

renders the soil more conductive to the effects of soil properties on infiltration processes, as discussed by Liu et al.⁹³. The third term, referred to as the snow severity index (SSI) in accordance with Diodato et al.⁹⁴, is instrumental in accounting for the contribution of snowmelt within the basin. This water input resulting from snowmelt serves as an effective supplier for GWR, demonstrating the capability to change the pattern of groundwater storage⁹⁵. The final term is the PDSI, which operates across both the current and preceding years.

**Extended model validation.** To validate the HHydroREM model in the absence of actual GWR data, we compared the model outputs with observed time-series data from the Boschetti well in Topino basin (Perugia, 43° 07' N, 12° 23' E, Umbria Region, for the upper TRB), along the period 1999 and 2020, and Peschiera discharge near Rieti (42° 24' N, 12° 42' E, Latium region, central TRB) between 1941 and 1974. In order to simulate the delay leading to groundwater recharge in both groundwater and spring discharge, we estimated the GWR(H) using a delayed drainage factor, GWR(H)ₖ over the preceding three years (y = -1, -2 and -3). This calculation was based on the simple moving average of the modelled discharge as:

$$GWR(H)_d = 0.3 \cdot GWR(H)_{y=-1} + 0.5 \cdot GWR(H)_{y=-2} + 0.3 \cdot GWR(H)_{y=-3} \quad (3)$$

In turn, the time-series of Eq. (3) was normalised. Then, in order to produce a hydrogeological response representative of the TRB, the Boschetti water table level and Peschiera discharge time-series were normalised and relatively averaged. The two time-series to be compared were normalised as:

$$\text{Stdz}_i = \frac{(\text{mean}(v_i) - v_i)}{\text{SD}_i} \quad (4)$$

where mean(vᵢ) is the mean value of time-series i, v is the value of the current year and SD is standard deviation of time-series i. In this way, the data were produced statistically and graphically, and the relative performance was provided.

**Model parameterisation.** The model's parameters were determined through a combination of statistical and physical evaluations, aiming for meaningful correlations. Estimating the value of the response variable was accomplished by discerning how changes in predictor data correspond to variations in the response variable, based on known values of the predictor variables. In the calibration sample, our initial objective was to gain an in-depth understanding of the components influencing the annual amount and temporal
distribution of groundwater recharge. Employing an iterative method and a trial-and-error approach, we identified the primary factors, providing a simplified explanation for the long-term dynamics of GWR. This involved a step-by-step process, with terms added or removed as necessary. Consequently, parameterisation underwent iterative adjustments to meet specific criteria of goodness-of-fit and error minimisation, as follows:

$$ \left\{ \begin{array}{l} R^2 = \max \\ |b - 1| = 0 \\ MAE = \min \end{array} \right. \qquad (5) $$

where the R² (optimum 1) is the goodness-of-fit of the linear regression of the calibration data to the estimates, $|b-1|$ is the difference from unity of the regression slope (b, optimum 1) between the calibration and modelled data, MAE (optimum 0) is the mean absolute error. The MAE extends the R² statistics to assess the accuracy of predictions in quantifying the amount of error. Skewness and kurtosis metrics were instead investigated to measure the data distribution and deviation from normality. To provide a comprehensive examination of the model's effectiveness, the two-sided Kolmogorov-Smirnov test was used to compare the distributions of Thorntwaite's GWR and predicted values, and the Durbin-Watson test was employed to investigate residual autocorrelation.

**Data analysis.** Data analysis was conducted using online statistical and graphical software tools, including STATGRAPHICS (http://www.statpoint.net/default.aspx), WESSA (https://www.wessa.net/tsa.wasp) and CurveExpert Professional 1.6 (https://www.curveexpert.net), to build and evaluate the model. Additionally, the Hurst (H) exponent (rate of chaos), indicating long-range dependence and cyclical-trend patterns in the time-series data, was estimated using SelQoS - self-Similarity and QoS analysis tool⁹⁶: 0.5 < H < 1.0 for positive memory (past trends tend to persist in the future); 0.0 < H < 0.5 for negative memory (past trends tend to reverting the future).

## Data availability

All data used in this study are freely available. Spatial patterns of global groundwater resources (Fig. 1) are freely available from https://www.whymap.org/whymap/EN/Home/whymap_node.html. The environmental, climatic and hydrogeological features of the study area (Fig. 2) are derived from various online sources, including: Mediterranean region and central Italian peninsula maps (Fig. 2a, b): https://www.clipsafari.com; mean annual precipitation distribution (Fig. 2c): https://opendata.dwd.de/climate/_environment/GPCC/html/fulldata-monthly_v2020_doi_download.html; Thorntwaite modified hydrological water balance (Fig. 2d): https://climatedataguide.ucar.edu/climate-data/gpcc-global-precipitation-climatology-centre and https://cran.uea.uea.ac.uk/cru/data/hrg/cru_ts_4.04; hydrological cycle dynamics and groundwater flow (Fig. 2e): https://www.noaa.gov/education/resource-collections/freshwater/water-cycle, https://iah.org/wp-content/uploads/2019/07/IAH_Climate-ChangeAdaptationGdwtr.pdf. The historical scenes of the Tiber valley and its river were sourced from the following online references (Fig. 3): https://www.kollerauktionen.ch/it/object-search.htm?term=FEODOR%20MICHALOWITSCH%20MATVEEV&origin= and https://it.wikipedia.org/wiki/Tevere#/media/File:The_Castel_Sant'Angelo_from_the_South.jpg. The standardised anomalies of the hydrogeological response (actual discharge) of the Tiber River Basin during 1941-1974 (Fig. 4d) can be accessed from http://www.bio.isprambiente.it/annalipdf. The data generated and analysed in this study are available at https://doi.org/10.5061/dryad.zs7h4jgs with a readme text file for use.

Received: 18 August 2023; Accepted: 18 January 2024; Published online: 29 January 2024

## References

1. Scanlon, B. R., Faunt, C. C., Longuevergne, L. & McMahon, P. B. Groundwater depletion and sustainability of irrigation in the US High Plains and Central Valley. *Proc. Natl. Acad. Sci. USA* **109**, 9320-9325 (2012).

2. Rodell, M. et al. Emerging trends in global freshwater availability. *Nature* **557**, 651–659 (2018).

3. Gampe, D. et al. Increasing impact of warm droughts on northern ecosystem productivity over recent decades. *Nat. Clim. Change* **11**, 772–779 (2021).

4. Pendergrass, A. G., Knutti, R., Lehner, F., Deser, C. & Sanderson, B. M. Precipitation variability increases in a warmer climate. *Sci. Rep.* **7**, 17966 (2017).

5. Persad, G. G., Swain, D. L., Kouba, C. & Ortiz-Partida, J. P. Inter-model agreement on projected shifts in California hydroclimate characteristics critical to water management. *Clim. Change* **162**, 1493–1513 (2020).

6. Fatichi, S., Peleg, N., Mastrotheodoros, T., Pappas, C. & Manoli, G. An ecohydrological journey of 4500 years reveals a stable but threatened precipitation-groundwater recharge relation around Jerusalem. *Sci. Adv.* **7**, eabe6303 (2021).

7. Letz, O. et al. The impact of geomorphology on groundwater recharge in a semi-arid mountainous area. *J. Hydrol.* **603**, 127029 (2021).

8. Long, D. et al. South-to-North water diversion stabilizing Beijing's groundwater levels. *Nat. Commun.* **11**, 3665 (2020).

9. Diodato, N., Bellocchi, G., Fiorillo, F. & Ventafredda, G. Case study for investigating groundwater and the future of mountain spring discharges in Southern Italy. *J. Mt. Sci.* **14**, 1791–1800 (2017).

10. Song, K., Ren, X., Mohamed, A. K., Liu, J. & Wang, F. Research on drinking-groundwater source safety management based on numerical simulation. *Sci. Rep.* **10**, 15481 (2020).

11. Diodato, N. et al. Climatic fingerprint of spring discharge depletion in the southern Italian Apennines from 1601 to 2020 CE. *Environ. Res. Commun.* **4**, 125011 (2022).

12. Stahl, K. et al. The challenges of hydrological drought definition, quantification and communication: an interdisciplinary perspective. *Proc. IAHS* **383**, 291–295 (2020).

13. Li, B., Rodell, M., Sheffield, J., Wood, E. & Sutanudjaja, E. Long-term, non-anthropogenic groundwater storage changes simulated by three global-scale hydrological models. *Sci. Rep.* **9**, 10746 (2019).

14. Fauzia, Surinaidu, L., Rahman, A. & Ahmed, S. Distributed groundwater recharge potentials assessment based on GIS model and its dynamics in the crystalline rocks of South India. *Sci. Rep.* **11**, 11772 (2021).

15. Jaafarzadeh, M. S., Tahmasebipour, N., Haghizadeh, A., Pourghasemi, H. R. & Rouhani, H. Groundwater recharge potential zonation using an ensemble of machine learning and bivariate statistical models. *Sci. Rep.* **11**, 5587 (2021).

16. Ouyang, Y., Zhang, J., Feng, G., Wan, Y. & Leininger, T. D. A century of precipitation trends in forest lands of the Lower Mississippi River Alluvial Valley. *Sci. Rep.* **10**, 12802 (2020).

17. Pfister, L., Savenije, H. H. G. & Fenicia F. Leonardo Da Vinci's Water Theory (IAHS Special Publication 9, 2009).

18. Immerzeel, W. W. et al. Importance and vulnerability of the world's water towers. *Nature* **577**, 364–369 (2020).

19. Jasechko, S. et al. The pronounced seasonality of global groundwater recharge. *Water Resour. Res.* **50**, 8845–8867 (2014).

20. Moeck, C. et al. A global-scale dataset of direct natural groundwater recharge rates: A review of variables, processes and relationships. *Sci. Total Environ.* **717**, 137042 (2020).

21. Lasagna, M., Mancini, S. & De Luca, D. A. Groundwater hydrodynamic behaviours based on water table levels to identify natural and anthropic controlling factors in the Piedmont Plain (Italy). *Sci. Total Environ.* **716**, 137051 (2020).

22. Pastore, N., Cherubini, C., Doglioni, A., Giasi, C. I. & Simeone, V. Modelling of the complex groundwater level dynamics during episodic rainfall events of a surficial aquifer in Southern Italy. *Water* **12**, 2916 (2020).

23. Behulu, F., Melesse, A. M. & Fiori, A. *Landscape Dynamics, Soils And Hydrological Processes In Varied Climates* (eds. Melesse, A. & Abtew, W.) 675–701 (Springer, 2016).

24. Caporali, E., Lompi, M., Pacetti, T., Chiarello, V. & Fatichi, S. A review of studies on observed precipitation trends in Italy. *Int. J. Climatol.* **41**, E1–E25 (2021).

25. Bencivenga, M. & Bersani, B. Influenza delle variazioni del clima sulle piene del Tevere a Roma. *Mem. Descr. Della Carta Geol. D’italia* **XCVI**, 377–385 (2014).

26. Maccioni, P., Kossida, M., Brocca, L. & Moramarco, T. Assessment of the drought hazard in the Tiber River Basin in central Italy and a comparison of new and commonly used meteorological indicators. *J. Hydrol. Eng.* **20**, 05014029 (2015).

27. Palmieri, S. et al. Hydrometeorological characterization of the Tiber Basin: role of evapotranspiration and soil storage in flood events. *Geol. Tec. Ambient.* **3**, 45–75 (2005).

28. Di Matteo, L., Dragoni, W., Maccari, D. & Piacentini, S. M. Climate change, water supply and environmental problems of headwaters: The paradigmatic case of the Tiber, Savio and Marecchia rivers (Central Italy). *Sci. Total Environ.* **598**, 733–748 (2017).

29. Durga Rao, K. H. V., Rao, V. V. & Dadhwal, V. K. Improvement to the Thorntwaite method to study the runoff at a basin scale using temporal remote sensing data. *Water Resour. Manag.* **28**, 1567–1578 (2014).
30. Mariotti, A. & Dell'Aquila, A. Decadal climate variability in the Mediterranean region: roles of large-scale forcings and regional processes. *Clim. Dyn.* **38**, 1129-1145 (2012).

31. Rohr, C., Camenisch, C. & Pribyl K. *European Middle Age In The Palgrave Handbook Of Climate History* (eds. White, S., Pfister C. & Mauelshagen, F.) 247-263 (Macmillan, 2018).

32. Shekhar, S. et al. Modelling water levels of northwestern India in response to improved irrigation use efficiency. *Sci. Rep.* **10**, 13452 (2020).

33. Sperna Weiland, F. C., van Beek, L. P. H., Kwadijk, J. C. J. & Bierkens, M. F. P. Global patterns of change in discharge regimes for 2100. *Hydrol. Earth Syst. Sci. Discuss.* **8**, 10973-11014 (2011).

34. Her, Y. et al. Uncertainty in hydrological analysis of climate change: multi-parameter vs. multi-GCM ensemble predictions. *Sci. Rep.* **9**, 4974 (2019).

35. Douville, H. & Willett, K. M. A drier than expected future, supported by near-surface relative humidity observations. *Sci. Adv.* **9**, eade625 (2023).

36. Qiu, J. et al. Nonlinear groundwater influence on biophysical indicators of ecosystem services. *Nat. Sustain.* **2**, 475–483 (2019).

37. Saccò, M. et al. Rainfall as a trigger of ecological cascade effects in an Australian groundwater ecosystem. *Sci. Rep.* **11**, 3694 (2021).

38. Zanchettin, D. et al. Atlantic origin of asynchronous European interdecadal hydroclimate variability. *Sci. Rep.* **9**, 10998 (2019).

39. Condon, L. E., Atchley, A. L. & Maxwell, R. M. Evapotranspiration depletes groundwater under warming over the contiguous United States. *Nat. Commun.* **11**, 873 (2020).

40. Manna, F., Walton, K. M., Cherry, J. A. & Parker, L. Five-century record of climate and groundwater recharge variability in southern California. *Sci. Rep.* **9**, 18215 (2019).

41. McCallum, J. L. et al. Assessing temporal changes in groundwater recharge using spatial variations in groundwater ages. *Water Resour. Res.* **56**, e2020WR027240 (2020).

42. Cuthbert, M. O. et al. Global patterns and dynamics of climate-groundwater interactions. *Nat. Clim. Change* **9**, 137-141 (2019a).

43. Diodato, N., Ljungqvist, C. F. & Bellocchi, G. A millennium-long climate history of erosive storms across the Tiber River Basin, Italy, from 725 to 2019 CE. *Sci. Rep.* **11**, 20518 (2021).

44. SIMN. Hydrological annals Servizio Idrografico and Mareografico Nazionale Italiano https://www.isprambiente.gov.it/it/progetti/cartella-progetti-in-corso/acque-interne-e-marino-costiere-1/progetti-conclusi/progetto-annali (1950-1975). (in Italian)

45. Civita, M. V. & Fiorucci, A. The recharge - discharge process of the Peschiera spring system (central Italy). *AQUA mundi* **02019**, 161–178 (2010).

46. Betocchi, A. *Del fiume Tevere* (Tipografia Elzeviriana, Rome, 1878).

47. Vaquero Piñeiro, M., 2001. Crescite incrociate: le piene del Tevere e lo sviluppo edilizio a Roma tra i secoli XVI e XVII in *I rischi del Tevere: modelli di comportamento del fiume di 10.1038/s41467-018-04207-7 Roma nella storia* (ed. Buonora, P.) 75-81 (Proceedings of the study seminar, 1998)

48. Tallaksen, L. M. & Van Lanen, H. A. *Hydrological Drought: Processes And Estimation Methods For Streamflow And Groundwater*. Volume 48 (Elsevier, 2004).

49. Hanel, M. et al. Revisiting the recent European droughts from a long-term perspective. *Sci. Rep.* **8**, 9499 (2018).

50. Markonis, Y., Hanel, M., Máca, P., Kyselý, J. & Cook, E. R. Persistent multi-scale fluctuations shift European hydroclimate to its millennial boundaries. *Nat. Commun.* **9**, 1767 (2018).

51. Büntgen, U. et al. Recent European drought extremes beyond Common Era background variability. *Nat. Geosci.* **14**, 190–196 (2021).

52. Lichtenthal, P. *Manuale di geografia fisica. Seconda Edizione Originale Rifusa E Considerevolmente Accresciuta* (Giuseppe Cioffi, 1854).

53. Cook, E. R. et al. Old World megadroughts and pluvials during the Common Era. *Sci. Adv.* **1**, e1500561 (2015).

54. Sadori, L. et al. Climate, environment and society in southern Italy during the last 2000 years. A review of the environmental, historical and archaeological evidence. *Quat. Sci. Rev.* **136**, 173-188 (2016).

55. Knight, J. R., Folland, C. K. & Scaife, A. A. Climate impacts of the Atlantic Multidecadal Oscillation. *Geophys. Res. Lett.* **33**, L17706 (2006).

56. Hurrell, J. W. Decadal trends in the North Atlantic Oscillation: regional temperatures and precipitation. *Science* **269**, 676–679 (1995).

57. Jones, P. D., Jónsson, T. & Wheeler, D. Extension to the North Atlantic Oscillation using early instrumental pressure observations from Gibraltar and South-West Iceland. *Int. J. Climatol.* **17**, 1433–1450 (1997).

58. Pinto, J. G. & Raible, C. C. Past and recent changes in the North Atlantic oscillation. *WIREs Clim. Change* **3**, 79–90 (2012).

59. Schoolman, E., Mensing, S. & Piovesan, G. From the late medieval to early modern in the Rieti basin (AD 1325–1601): paleoecological and historical approaches to a landscape in transition. *Hist. Geogr.* **46**, 103–128 (2018).

60. Polovodova Asteman, I., Filipson, H. L. & Nordbäran, K. Tracing winter temperatures over the last two millennia using a north-east Atlantic coastal record. *Clim. Past* **14**, 1097–1118 (2018).

61. Borghi, A. *DEscrizione Geografica, Fisica E Naturale Del Lago Trasimeno Comunemente Detto Lago Di Perugia* (Tipografia Bassoni, 1821).

62. Gasperini, L. et al. Late Glacial and Holocene environmental variability, Lago Trasimeno, Italy. *Quat. Int.* **622**, 21–35 (2022).

63. Eddy, J. A. The maunder minimum. *Science* **192**, 1189–1202 (1976).

64. Frosini, P. Il lago Trasimeno e il suo antico emissario. *Boll. Soc. Geogr. It., Serie VIII* **11**, 6–15 (1958). (in Italian).

65. De Sanctis, G. *Dizionario Statistico De' Paesi Del Regno Delle Due Sicilie* (Nabu Press, 1840).

66. Tozzi, I. *Il Culto Delle Madonne Arboree Nel Territorio Delle Diocesi di Rieti e Sabina* (Ileana, 2003).

67. Ionita, M., Dima, M., Nagavciuc, V., Scholz, P. & Lohmann, G. Past megadroughts in central Europe were longer, more severe and less warm than modern droughts. *Commun. Earth Environ.* **2**, 61 (2021).

68. Burzigotti, R., Dragoni, W., Evangelisti, C. & Gervasi, L. *The Role Of Lake Trasimeno (Central Italy) In The History Of Hydrology And Water Management* (International Water History Association, 2003).

69. García-Ruiz, J. M., López-Moreno, J. I., Vicente-Serrano, S. M., Lasanta-Martinez, T. & Beguería, S. Mediterranean water resources in a global change scenario. *Earth-Sci. Rev.* **105**, 121–139 (2011).

70. Yin, X.-A., Yang, X. & Yang, Z.-F. Using the R/S method to determine the periodicity of time series. *Chaos Solitons Fractals* **39**, 731–745 (2009).

71. Sheng, H. & Chen, Y. Q. *Robustness Analysis Of The Estimators For Noisy Long-range Dependent Time Series* (American Society of Mechanical Engineers, 2009).

72. Quian, B. & Rasheed, K. Hurst exponent and financial market predictability. in *Proceedings of the 2nd IASTED International Conference on Financial Engineering and Applications* (ed. Hamza, M. H.) 203-209 (International Association of Science and Technology for Development, 2004).

73. Hale, G. E., Ellerman, F., Nicholson, S. B. & Joy, A. H. The magnetic polarity of sunspots. *Astrophys. J.* **49**, 153–178 (1919).

74. Zanchettin, D., Rubino, A., Traverso, P. & Tomasino, M. Impact of variations in solar activity on hydrological decadal patterns in northern Italy. *J. Geophys. Res.* **113**, D12102 (2008).

75. Ojala, A. E. K., Launonen, I., Holmström, L. & Tiljander, M. Effects of solar forcing and North Atlantic oscillation on the climate of continental Scandinavia during the Holocene. *Quat. Sci. Rev.* **112**, 153–171 (2015).

76. ISTAT. Statistics on water resources, water use and wastewater treatment – Development of data collections systems and statistical methods for indicators at the sub-national level. *Social and Environmental Statistics Department Socio-demographic and Environmental Statistics Directorate State of Environment Division Water Resources And Climate Unit https://ircabc.europa.eu/sd/a/bd29a96c-09ff-4668-864d-d02685f7767/1stat_grant_final_report_water.pdf* (2010).

77. *Şen, Z. Applied Drought Modeling, Prediction, And Mitigation* (ed. Şen, Z.) 321-391 (Elsevier, 2015).

78. Kommadath, A. Estimation of natural ground water recharge. *Ground Water Hydrogeol.* **51**, 80–91 (2000).

79. Scanlon, B. R., Healy, R. W. & Cook, P. G. Choosing appropriate techniques for quantifying groundwater recharge. *Hydrogeol. J.* **10**, 18–39 (2002).

80. Thornthwaite, C. W. An approach toward a rational classification of climate. *Geogr. Rev.* **38**, 55–94 (1948).

81. Mather, J. R. *Use Of The Climatic Water Budget To Estimate Streamflow In Use Of The Climatic Water Budget In Selected Environmental Water Problems* (ed. Mather, J. R.) **32**, 1–52 (Elmer, N. J., C.W. Thornthwaite Associates, 1979).

82. McCabe, G. J. & Wolock, D. M. Future snowpack conditions in the western United States derived from general circulation model climate simulations. *J. Am. Water Resour. Assoc* **35**, 1473–1484 (1999).

83. Vandewiele, G. L. & Win, N. L. Monthly water balance models for 55 basins in 10 countries. *Hydrol. Sci. J.* **43**, 687–699 (1998).

84. Mammoliti, E., Fronzi, D., Mancini, A., Valigi, D. & Tazioli, A. WaterbalANce, a WebApp for Thornthwaite-Mather water balance computation: comparison of applications in two European watersheds. *Hydrology* **8**, 34 (2021).

85. McCabe, G. J. & Markstrom, S. L. A monthly water-balance model driven by a graphical user interface. *Open-File report 2007-1088* (U.S. Geological Survey, 2007).

86. Landerer, F. W. & Swenson, S. C. Accuracy of scaled GRACE terrestrial water storage estimates. *Water Resour. Res.* **48**, W04531 (2012).

87. Van Loon, A. F. Hydrological drought explained. *WIREs Water* **2**, 359–392 (2015).

88. Oke, I. A. Reliability and statistical assessment of methods for pipe network analysis. *Environ. Eng. Sci.* **24**, 1481–1489 (2007).

89. Cuthbert, M. O. et al. Observed controls on resilience of groundwater to climate variability in sub-Saharan Africa. *Nature* **572**, 230–234 (2019b).

90. Diodato, N. & Fiorillo, F. Complexity-reduced in the hydroclimatalogical modelling of aquifer's discharge. *Water Environ. J.* **77**, 170–176 (2012).

91. Esit, M. et al. Seasonal to multi-year soil moisture drought forecasting. *npj Clim. Atmos. Sci.* **4**, 16 (2021).
92. Goni, I. B. et al. Groundwater recharge from heavy rainfall in the southwestern Lake Chad Basin: evidence from isotopic observations. *Hydrol. Sci. J.* **66**, 1359–1371 (2021).

93. Liu, Y., Cui, Z., Huang, Z., López-Vicente, M. & Wu, G.-L. Influence of soil moisture and plant roots on the soil infiltration capacity at different stages in arid grasslands of China. *Catena* **182**, 104147 (2019).

94. Diodato, N., Büntgen, U. & Bellocchi, G. Mediterranean winter snowfall variability over the past Millennium. *Int. J. Climatol.* **39**, 384–394 (2019).

95. Wu, W. Y. et al. Divergent effects of climate change on future groundwater availability in key mid-latitude aquifers. *Nat. Commun.* **11**, 3710 (2020).

96. Ramírez-Pacheco, J., Torres-Román, D., Toral-Cruz, H. & Vargas, L. *High-performance Tool For The Test Of Long-memory And Self-similarity In Simulation Technologies In Networking And Communications* (eds. Khan Pathan, A.-S., Monowar, M. M. & Khan, S.) 93-114 (CRC Press, 2015).

97. Banzato, F. et al. Relationship between rainfall and water table in a coastal aquifer: the case study of Castelporziano presidential estate. *Ital. J. Groundw.* **8**, 27–33 (2019).

98. Giraudi, C. Late Pleistocene and Holocene lake-level variations in Fucino Lake (Abruzzo, central Italy) inferred from geological, archaeological and historical data. *Palaoklimaforschung* **25**, 1–17 (1998).

99. Mann, M. E. et al. Global signatures and dynamical origins of the Little Ice Age and Medieval Climate Anomaly. *Science* **326**, 1256–1260 (2009).

100. PAGES 2k Consortium. Continental-scale temperature variability during the past two millennia. *Nat. Geosci.* **6**, 339–346 (2013).

## Acknowledgements

Nazzareno Diodato and Gianni Bellocchi performed this research as an investigator-driven study without financial support.

## Author contributions

Nazzareno Diodato conceptualised and developed the original research design, and collected and analysed the historical documentary data. Gianni Bellocchi co-wrote the article and contributed to the interpretation of findings. All authors critically reviewed and approved the final manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-024-01229-6.

**Correspondence** and requests for materials should be addressed to Gianni Bellocchi.

**Peer review information** *Communications Earth & Environment* thanks Alejandro Fernandez and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Patricia Spellman and Aliénor Lavergne. A peer review file is available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024