# How a millennium-long interplay of controlling climate factors over Southern Europe (1054–2020 CE) has driven the downward trend in Eastern Po Plain water reserves

Nazzareno Diodato, Fausto Tomei, Gianni Bellocchi

► To cite this version:

Nazzareno Diodato, Fausto Tomei, Gianni Bellocchi. How a millennium-long interplay of controlling climate factors over Southern Europe (1054–2020 CE) has driven the downward trend in Eastern Po Plain water reserves. Water International, 2025, pp.1-23. 10.1080/02508060.2025.2497216. hal-05100504

HAL Id: hal-05100504

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

Submitted on 6 Jun 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
# How a millennium-long interplay of controlling climate factors over Southern Europe (1054–2020 CE) has driven the downward trend in Eastern Po Plain water reserves

Nazzareno Diodato, Fausto Tomei and Gianni Bellocchi

## QUERY SHEET

This page lists questions we have about your paper. The numbers displayed at left are hyperlinked to the location of the query in the paper.

The title and author names are listed here as they will appear in your paper and the Table of Contents. Please check that this information is correct and let us know if any changes are needed. Also check that affiliations, funding information and conflict of interest statements are correct.

Please review your paper as a whole for typos and essential corrections. Note that we cannot accept a revised manuscript at this stage of production or major corrections, which would require Editorial review and delay publication.

## AUTHOR QUERIES

**Q1** Please provide missing department for the affiliation c.

**Q2** Please note that the year 2016 in the mention (citation) of reference "Zanchettin et al., 2016" has been changed to 2019 to match the reference list. Check and if necessary, correct directly in the text or reference.

**Q3** Please give "SLP" in full.

**Q4** The disclosure statement has been inserted. Please correct if this is inaccurate.

**Q5** The PubMed (http://www.ncbi.nlm.nih.gov/pubmed) and CrossRef (www crossref. org/) databases have been used to validate the references. Mismatches between the original manuscript and PubMed or CrossRef are tracked in red font. Please provide a revision if the change is incorrect. Do not comment on correct changes.

**Q6** Please provide missing issue number and DOI or URL, if available, for the "Conte, 1992" references list entry.

**Q7** Please provide missing page number for the "Finzi, 1986" references list entry.

**Q8** The Reference "He, 2021" is listed in the references list but is not cited in the text. Please either cite the reference or remove it from the references list.

**Q9** Please provide missing issue number and DOI or URL, if available, for the "Kumar, 2012" references list entry.

**Q10** Please provide missing issue number and DOI or URL, if available, for the "Pilli et al., 2012" references list entry.
Q11 The Reference "Yang, 2016" is listed in the references list but is not cited in the text. Please either cite the reference or remove it from the references list.

Q12 Please provide the high-resolution quality for All figures.

Q13 The year for "Thornthwaite, 1941" has been changed to 1948 to match the entry in the references list. Please provide revisions if this is incorrect.

Q14 Please check whether the affliaton has been set correctly.
RESEARCH ARTICLE

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

Check for updates

# How a millennium-long interplay of controlling climate factors over Southern Europe (1054–2020 CE) has driven the downward trend in Eastern Po Plain water reserves

Nazzareno Diodato<sup>a</sup>, Fausto Tomei<sup>b</sup> and Gianni Bellocchi<sup>a,c</sup>

Q14
Q1

## ABSTRACT

Groundwater can unveil climate signals and improve long-term reconstructions, but historical studies often overlook groundwater depth due to lack of data. This study presents the longest annual groundwater depth time series in southern Europe (1054–2020 CE) from a relatively deep aquifer in the Eastern Po Plain, northern Italy. The record links groundwater depth oscillations around its mean value (–147 cm ±19 standard deviation) to Atlantic climate patterns and water balance. A significant downward trend after about 1606 indicates a depletion of regional water reserves and underlines the need to strengthen groundwater management to address potentially worsening drought trends in southern Europe.

## ARTICLE HISTORY

Received 2 April 2024

Accepted 21 April 2025

## KEYWORDS

Climate signals;

groundwater; long-term

reconstruction; Southern

Europe

5

10

15

20

## Introduction

The Sustainable Development Goals (SDGs), adopted by the United Nations in 2015 to achieve global sustainability by 2030 (https://sdgs.un.org/2030agenda), recognize the critical role of groundwater in successfully achieving multiple goals. Groundwater's significance is evident across various SDGs, including promoting sustainable groundwater use in agriculture (SDG 2), safeguarding human health and preventing waterborne diseases (SDG 3), ensuring clean water and sanitation (SDG 6), emphasizing the importance of groundwater in urban water supply and development (SDG 11), recognizing groundwater's contribution to climate resilience (SDG 13) and understanding groundwater's impact on marine ecosystems (SDG 14).

25

30

While nowadays there is a growing attention to the significance of groundwater resources, it is essential to acknowledge that its importance has been recognized since ancient times. Groundwater has played a crucial role in sustaining various ecosystems and meeting human needs throughout history (Rodell et al., 2018).

Historical figures like Bernard Palissy (ca. 1510–ca. 1589), French artist and naturalist, in analogy with what the Italian High Renaissance polymath Leonardo da Vinci

35

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

<footnote><sup>a</sup>Met European Research Observatory, International Affiliates Program of the University Corporation for Atmospheric Research, Benevento, Italy; <sup>b</sup>Osservatorio clima, Arpae Emilia-Romagna, Bologna, Italy;</footnote>

<footnote><sup>c</sup>Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France</footnote>
<figure><img src="image_4.png" /></figure>

(1452–1519) had sketched (De Lorenzo, 1939), had already identified its importance in his time (Morley, 1853, p. 324, Supplementary file).

From historical figures like Palissy and da Vinci to contemporary concerns around the SDGs, the importance of groundwater has been recognized for centuries. As we face global challenges, it is imperative that we understand and protect this vital resource for future generations. However, despite their significance, studies on subsurface hydrology are still primarily limited to short timescales (Li et al., 2019). As Endo (2022) pointed out, groundwater management has received much less attention than surface-water management in the fields of environmental and water history, as obtaining homogeneous groundwater data poses challenges. Information from monitoring wells is sometimes intermittent, irregular, infrequent or limited to brief periods (Evans et al., 2020), hindering investigations on historical groundwater dynamics (Fauzia et al., 2021). Overcoming this obstacle will be crucial in better understanding and effectively managing this critical resource.

Figure 1 illustrates the complexity of hydrogeological structures globally and their implications for groundwater recharge. In Southern Europe, this complexity is particularly pronounced, with aquifers exhibiting intricate formations that require careful planning and management due to the limited information available on these resources

40

45

50

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Global groundwater recharge and main hydrogeological structures. This map provides information on the quantity, quality and vulnerability of Earth's groundwater resources, and helps experts, decision-makers and the public to understand, address and communicate groundwater-related issues. The southern Europe region (highlighted in the red rectangle) is supplied by aquifers with complex hydrogeological structures, requiring careful planning and management due to limited information available on this hidden resource. The map is sourced from the World-wide Hydrogeological Mapping and Assessment Programme (WHYMAP) of the German Federal Institute for Geosciences and Natural Resources and the United Nations Educational, Scientific and Cultural Organization (https://www.whymap.org/whymap/EN/Home/whymap_node.html).</figcaption></figure>
<figure><img src="image_6.png" /></figure>

(Cramer et al., 2018). Understanding these hydrogeological features is crucial to effec-
tively address groundwater recharge issues, especially in regions where water scarcity
may occur.

55

Groundwater storage (size measured by its depth) is the precipitation remain-
ing after canopy interception, minus the sum of evapotranspiration and runoff,
which are affected by change and variability in climate, vegetation and soil
properties (Kim & Jackson, 2012). Temperature influences the hydrologic cycle
by increasing evaporation of usable surface water and plant transpiration. In turn,
these changes can impact precipitation volumes, timing and intensity rates, as well
as surface and groundwater flows and storage, which are essential for the survival
of streams, lakes, wetlands and other ecosystems that depend on them, from
regional to local scales.

60

65

Variations in these factors control the amount of water that reaches the surface, evaporates or seeps into the atmosphere, accumulates as snow or ice, infiltrates into the groundwater system, flows out of the ground and finally becomes the baseflow of streams and rivers (Kumar, 2012). In this way, aquifer recharge in Italy is particularly sensitive to climate change, putting them in a critical situation in case of prolonged drought (Di Matteo et al., 2021) and with increasing rainfall extremes (Treppiedi et al., 2021). In Italy groundwater reserves are mostly unprotected from the consequences of climate change brought about by climate variability combined with increased water demand.

70

75

In this context, historical reconstruction of groundwater resource data remains a challenge for environmental hydrogeology (McCallum et al., 2020). Ultimately, **because** predictions remain model- and scenario-dependent (Wu et al., 2020), global hydrology models are ineffective for planning of water management strategies, which can only be achieved and implemented locally at the aquifer scale. Also past hydroclimatic changes and, in particular, groundwater-climate interactions are poorly understood at this scale (Saccò et al., 2021), not least because of the difficulty of reconstructing groundwater data over long periods (Zanchettin et al., 2019).

80

Q2

A

To better invest in the future, we need to improve our understanding of how ground-
water depth has responded to historical climatic changes, beyond human impacts (Li
et al., 2019). This study builds on the concept and methods used in a simple groundwater
depth model to bridge the gap between the large amounts of knowledge exchanged in
climate research and the paucity of data for historical modelling. For this purpose, the
Groundwater REsource Model (GREM) was developed to estimate past values of annual
groundwater depth for the Eastern Po Plain, where groundwater depth is supported by
moderate amounts of precipitation. The proposed model is sober because it relies on
simplified inputs including large-scale teleconnection indices, a regional water balance
proxy and a subregional storm intensity index. These indicators were chosen to over-
come the limitations imposed by more complex, physics-based models, which would not
have allowed reconstructing the groundwater depth time series in historical times due to
the lack of past hydrogeological input data. In this way, GREM helps to better understand
how groundwater interacts with the regional water cycle and responds to the complex
vicissitudes of how past climate sustained the Eastern Po Plain aquifers.

85

90

95

Here, we addressed this critical issue, thereby contributing to research on groundwater
and climate change in the Eastern Po Plain, northern Italy (Figure 2(a)).

100
<figure><img src="image_7.png" /></figure>

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

**Figure 2.** Environmental setting and hydrological response in the study area. (a) Mean annual precipitation in Italy from 1971 to 2000, sourced from Esposito et al. (2015) and provided by the Italian National Institute for Environmental Protection and Research (https://www.isprambiente.gov.it/pre_meteo/idro/BIGBANG_CC_ISPRA.html). (b) Mean annual potential groundwater recharge in northern Italy from 1981 to 2000, with specific locations indicated by dots. The Cadriano piezometer, used for calibration, is denoted by a small black dot. Additionally, the T12-Tonale (source: http://www.comune.giussano mb.it/flex/cm/pages/ServeAttachment.php/L/IT/D/9%252F7%252F6%252FD.c78fd0d1483fb4568af8/P/BLOB%3AID%3D1500) and T16-Tortona piezometers, used for extended evaluation, are shown in red dots. These data are courtesy of Professor Simon Stisen from the Department of Hydrology at the Geological Survey of Denmark and Greenland. (c) The modified Thornthwaite (1948) hydrological water balance over the Po Plain for the period 1981–2010. This analysis is based on an updated version of the WebWIMP programme – WebWIMP 1.01 – as implemented by the University of Delaware, USA, in 2003 (http://climate.geog.udel.edu/~wimp/index.html).

Q13

# Environmental setting and methods

## Geographical and climatic characteristics of the study area

With a length of 652 km, the Po River is the longest river in Italy. The Po Plain (also known as Po Valley, Plain of the Po or Padan Plain) is a major geographical feature of Northern Italy around the course of the Po River that projects into the Adriatic Sea. The Po Plain extends in a west-east direction from the Western Alps to the Adriatic Sea, with an area of ~46,000 km² including its extension in the Venetian-Friulian Plain. The latter
<figure><img src="image_9.png" /></figure>

is not actually related to the Po basin, as its flatlands do not drain into the Po, but effectively combine into an uninterrupted plain, making it the largest in southern Europe. The Eastern Po Plain territory (our study area) falls within the wide alluvial plain of the lower basin of the Po River.

110

The area is located in the altitude range from 0 to 100 m a.s.l., with gentle slopes. The typical climate characteristics of the Po Plain in the northern part and those of the Apennine chain, in the southernmost part, affect the area. While, in general, the groundwater recharge potential in the Po Plain is relatively low (Figure 2(b)), in winter and early spring the groundwater recharge is sustained in the Eastern Po Plain (Figure 2(c)). At the end of the spring and during the summer months, rainfall is partly associated with convective phenomena and characterized by intense spatial variability and greater water deficit. Autumn is the season when water reserves are replenished before recharging becomes sustained again in early winter (Figure 2(c), blue). Together with climate change, human interactions have recently modified the water cycle by draining groundwater and surface water, and modifying the land cover to meet the demand in the Po Plain.

115

120

For modelling purposes, the calibration time series was obtained from the Department of Agricultural and Food Sciences of the University of Bologna and University of Bologna Experimental Farm for the period 1975-2012, and referred to the piezometric station of Cadriano (44.55 N, 11.40 E, 27 m a.s.l.; Figure 2(c)), in the metropolitan area of Bologna.

125

The Cadriano area collects ~740 mm yr⁻¹ (mean for the period 1971-2010), with the surrounding areas exceeding 1000 mm yr⁻¹ (Figure 2(a)). The site was chosen to develop and parameterize the GREM because of the availability here of annual piezometric station measurements (groundwater depth) from 1975 to 2012 and because it is deep enough to prevent the annual time series of groundwater depth data from being spoilt by the influence of human activities. The reliability of the GREM estimates was corroborated by an evaluation extended to the Western Po Plain (Figure 2(c)). For the piezometric station T16-Tortona (44.46 N, 8.86 E., 230 m a.s.l.), near the city of Alessandria, the time series was derived from Lasagna et al. (2019) for the period 2002-2017. For the piezometric station T12-Tonale (45.70 N, 9.20 E, 266 m a.s.l), in the municipality of Giussano (near Milan), the time series was obtained from Comune di Giussano for the period 1977-1998 (http://www.comune.giussano mb.it/flex/cm/pages/ServeAttachment.php/L/IT/D/9%252F7%252F6%252FD.c78fd0d1483fb4568af8/P/BLOB%3AID%3D1500). In this way, the modelling study was designed to calibrate the model on the Cadriano aquifer, for which a large data set was available, and then to evaluate it on two independent data sets of different type: the shallow unconfined aquifer of T16-Tortona and the deep, partly confined aquifer of T12-Tonale. While one of the independent data sets came from a similar well (Tortona), the second (Tonale) was from a deep aquifer (in which the observed variable is referred to the piezometric level in m a.s.l.), which is little vulnerable to contamination and suitable for drinking. As a paleochannel of the Lambro River (a left tributary of the Po, 130 km long), it is a depression filled almost exclusively with gravelly and conglomeratic sediments, bounded to the east by the rocky substrate and to the west by clays in marine facies, forming a monolayer aquifer. The Lambro River has its source below Lake Como, north of the city of Erba (45.82 N, 09.22 E), and has been diverted to Lake Pusiano (45.80 N, 9.27 E), to avoid flooding the villages downstream (https://web.archive.org/web/20160712175908/http://www.cavodiotti.it/cenni.html). The

130

135

140

145

150
<figure><img src="image_10.png" /></figure>

Lambro paleochannel, bounded by the rocks on the sides (east and west), produces a large and deep aquifer, fed by upstream recharges upstream. Although the evaluation of the model in this well may be critical, we maintained it in the study to test the model's potential in capturing situations with longer recharge times (up to months and years) and outflows that are predominantly anthropogenic, but mainly related to climatic conditions (the warmer it is, the higher the water consumption). In fact, the first measurements, relative to the Tonale monitoring network, show a rise in the water table due to the intense above-average rainfall that occurred from the end of 1976 to the first part of 1978, with the water table reaching its minimum distance from the ground level in the summer of 1978. After 1979, there was a gradual decrease, with positive fluctuations due to the considerable precipitation, including snowfall, in the winter of 1985 and various autumn periods, particularly those with heavy rainfall in 1993 and 1996, the effects of which were felt until the first months of the following years (relative highs in February 1994 and February 1997). Once parameterized and evaluated, the model was applied to estimate the long-term groundwater depth of the Eastern Po Plain through distinct climatic periods, here defined as the Medieval Climate Anomaly (MCA, 801–1299 CE), the Little Ice Age (LIA, 1300–1849 CE) and the Modern Warming Period (MWP, 1850–2020 CE). In doing so, we have produced and discussed a long time series of annual groundwater depth data, unique in its length (1054–2020 CE) and thus valuable for hydrological studies in southern Europe.

The North Atlantic Oscillation and Atlantic Multidecadal Variability time series, used as model inputs, were derived from Cook et al. (2015) and Wang et al. (2017), respectively. The North Atlantic Oscillation is the most important pattern influencing weather and climate variability and precipitation extremes in Europe, from synoptic to millennial timescales (Casanueva et al., 2014). It consists of a mean sea-level pressure dipole pattern between the subtropical and subpolar areas of the North Atlantic and characterizes the intensity and location of the westerly flow. On multidecadal timescales, the Atlantic Multidecadal Variability, which is a warming/cooling pattern of sea-surface temperatures in the North Atlantic with periodicities of ~60 years, is also known to influence European surface temperatures and regimes over Europe, including northern Italy (Sutton & Dong, 2012). Other inputs are the Annual Storm Severity Index from Diodato et al. (2020) and the Palmer Drought Severity Index from Cook et al. (2015). The Annual Storm Severity Index Sum exemplifies storm anomalies as an annual sum of monthly integer scores from 1 (stormy event) to 4 (extraordinary storm event). The Palmer Drought Severity Index is a standardized regional index, not explicitly bounded, which is based on a simplified soil water balance and estimates relative soil moisture conditions, typically ranging from – 4 (extreme drought) to +4 (extreme wetness).

## Groundwater Resource Model

The historical estimation of the annual mean groundwater depth (cm) for the Eastern Po Plain was performed using a non-linear multivariate regression model, hereafter referred as GREM. This model integrates several key climatic and hydrological variables, each selected for their observed correlation with groundwater variations. The inclusion of these variables aligns with established hydrological principles and extensive observations, and supports the model's alignment with regional groundwater behaviour (Beven, 2012).
<figure><img src="image_11.png" /></figure>

Extensive research has shown that groundwater oscillations are correlated with climatic drivers, such as teleconnection indices (Rust et al., 2022) and hydrological forcing (Singh & Borro, 2019). Consequently, we have used a parsimonious approach to derive a non-linear equation that incorporates some of these influential drivers.

The model, introduced in Equation (1), is formulated as a function $-Z = A(X + Y + Z)$ – of three primary variables (X, Y and Z), representing different spatial scales and hydrological drivers: large-scale teleconnective climatic factor (LSTCF), a regional drought function (RSDF) and a subregional storm severity factor (SRSHF). This approach enables the model to capture both broad climatic patterns and localized hydrological responses. The equation is as follows:

$$GREM_{(y=0)} = A \cdot (f(LSTCF) + f(RSDF) + f(SRSHF)) + B \quad (1)$$

where GREM$_{(y=0)}$ is the estimated annual mean groundwater depth (cm) for the current year (y = 0), and A and B are constants derived through regression fitting. Explaining the three main terms, the following equation is obtained for GREM$_{(y = 0)}$:

$$GREM_{(y=0)} = A \cdot \left[ \text{Prc}_{25}(\text{NAO} + \text{AMV})^{y-5}_{y-1} + \text{PDSIs} + \left( 1 + \sum_{y-0}^{y-2} \text{ASSIS} \right)^{\prime\prime} \right] + B \quad (2)$$

where NAO and AMV are the North Atlantic Oscillation and Atlantic Multidecadal Variability, respectively; ASSIS is the Annual Storm Severity Index Sum and PDSI is the Palmer Drought Severity Index.

The following terms clarify each variable and its physical significance:

(1) Large-scale teleconnective climatic factor (LSTCF). This term combines the North Atlantic Oscillation and Atlantic Multidecadal Variability indices, two climate oscillations that significantly affect weather and precipitation patterns across southern Europe (De Vita et al., 2022; Rust et al., 2012). Deeper aquifers respond slowly to large-scale climate shifts but are more sensitive to short-term climate variability (Clifton et al., 2010). Conversely, shallow groundwater systems, like those in the Eastern Po Plain's unconsolidated sediments or fractured bedrock, respond rapidly to smaller-scale climate variability (Kundzewicz & Döll, 2009). In particular, a negative North Atlantic Oscillation phase brings mild, rainy conditions associated with Atlantic weather patterns (advective precipitation), while a negative Atlantic Multidecadal Variability phase results in rainy intervals due to Mediterranean cyclones sustained by cold, unstable atmospheric conditions (Conte, 1992). For this model, we use the 25th percentile of North Atlantic Oscillation and Atlantic Multidecadal Variability coupling over a five-year period prior to the current year of estimation, within the range of -0.5 to +0.5 (Figure 3).

This aligns with the findings of Neves et al. (2019) in Portugal, who observed the strongest covariation occurring in the 6–10 years band in the response of North Atlantic Oscillation to piezometric levels. Within the Easter Po Plain region, Meli and Romagnoli (2022) also identified partial correlations between the North Atlantic Oscillation and Atlantic Multidecadal Oscillation with rainfall and river discharge. In addition, the percentile (Prc) correlation of the groundwater depth not only improved
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 3.** Schematic illustration of groundwater depth variation with percentile analysis. The figure demonstrates the variation in groundwater depth represented by the 25th percentile (Prc25) of the combined North Atlantic Oscillation (NAO) and Atlantic Multidecadal Variability (AMV) oscillations over the five-year window prior to the current year ($y = 0$). This relationship emerged using the percentile statistical operator, distinct from other operators.</figcaption></figure>

model performance but also led to a new interpretation of the relationship between the behaviour of the Po aquifer and the combination of the climatic teleconnective indices. Notably, groundwater exhibits a positive response, especially when the coupling of the indices is negative, but also at times when it is positive, indicating a non-linear teleconnective relationship (Figure 3). The only neutral coupling of the indices, with values around zero, corresponds to some of the lowest groundwater depth values, approximately –150 m.

240

(2) *Regional-scale drought function (RSDF).* This term incorporates the Palmer Drought Severity Index, a widely used metric that reflects soil water budget and provides a broad assessment of regional conditions. Positive Palmer Drought Severity Index values suggest adequate soil moisture, which supports recharge, whereas negative Palmer Drought Severity Index values indicate drought-induced evapotranspiration, reducing soil and groundwater reserves. The Palmer Drought Severity Index term is essential because drought impacts extend across water supply, agriculture, socioeconomic development and ecosystem health (Pokhrel et al., 2021), making it a critical factor for sustainable groundwater management (Hanel et al., 2018). These conventional indices have been used in meteorological, agricultural and hydrological droughts and in drought monitoring and projection (Garamhegyi et al., 2018). This could mean that there is no or low groundwater deficit with positive Palmer Drought Severity Index values. Conversely, in warm-season drought events, a rapid increase in excessive evapotranspiration is observed, which will deplete soil moisture and, in turn, groundwater depth, with negative Palmer Drought Severity Index values.

245

250

255

(3) *Subregional storm severity hydrological factor (SRSHF).* The SRSHF term, based on the Annual Storm Severity Index Sum from Diodato et al. (2020), which takes into account both the duration and intensity of precipitation, the amounts of precipitation fallen in year $y = 0$ (current year), and in two previous years $y = -1$ and $y = -2$, identifying the memory effect that accumulates in the current year $y = 0$. This term

260
is a modified empirical relationship from Oke (2007), between the groundwater recharge and the corresponding monsoon storm values using the non-linear regression equation. This multiyear memory is particularly relevant in the Eastern Po Plain, where storm-driven recharge can mitigate deficits even in years with lower overall precipitation (Cuthbert et al., 2019; Esit et al., 2021).

In summary, we selected each term in our model not only for its statistical significance but also for its physical relevance to groundwater recharge processes. The model structure reflects the interplay between seasonal climate patterns and local hydrological events. Broad-scale climate indices, such as the North Atlantic Oscillation and Atlantic Multidecadal Variability, capture the influence of winter precipitation, a key driver of groundwater recharge in temperate regions. The Palmer Drought Severity Index integrates soil moisture conditions over larger areas, while the Annual Storm Severity Index Sum provides additional detail on how extreme precipitation events can directly recharge aquifers, even during dry periods. This parsimonious, non-linear regression approach provides a meaningful representation of groundwater dynamics across different temporal and spatial scales. By accounting for both large-scale climate teleconnections and local hydrological variability, the GREM model provides a robust, data-driven framework for estimating groundwater levels and supporting informed water management decisions.

## Model parameterization and data analysis

The parameterization of the model was inferred through a combination of statistical and physical evaluations, incorporating observed correlations. This allowed the model to be ready for prediction run. To achieve this, we initially employed a lumped model approach, using trial and error as the most basic optimization method. During this process, the model parameters were adjusted, and the fit between the model results and the calibration data set was measured. This iterative process was repeated multiple times to ensure that the observed and predicted values converged as closely as possible to provide reliable predictions. The modelling routine underwent several iterations to meet the following conditions:

$$R^2 = \max; |b-1| = \min; \text{MAE} = \min \quad (3)$$

R² (optimum 1) represents the goodness-of-fit of the linear regression of the actual data to the estimates. Its optimum value is 1, indicating a perfect fit between the observed data and the model estimates. On the other hand, |b-1| represents the difference between the actual and modelled data regression slope (b). The optimum value for this parameter is 1, which indicates that the model's regression slope is exactly the same as the actual data. Lastly, MAE (mean absolute error) is a metric used to quantify the accuracy of the model's predictions. Its optimum value is 0, which means that the model's predictions exactly match the observed data with no error. These metrics are important indicators of the model's performance and how well it aligns with the actual data. Higher values of R² and |b-1| closer to 1, along with a lower MAE closer to 0, indicate a better-performing model with higher accuracy and goodness-of-fit.
<figure><img src="image_14.png" /></figure>

Spreadsheet-based exploratory data analysis and model building was performed with the help of the online statistical software STATGRAPHICS (http://www.statpoint.net/default.aspx), with graphical support from AnClim (http://www.climahom.eu/software-solution/anclim; Štěpánek, 2005) and CurveExpert Professional 1.6 (https://www.curveexpert.net).

305

# Results and discussion

## Model parameterization and evaluation

Groundwater depth data obtained from measurements at the Cadriano piezometer station for the period from 1975 to 2008 was used as a calibration data set. We discarded two values in the years 1989 (–185 cm) and 2002 (–170 cm), which were below the lower 95% confidence limit and were considered outliers (orange dots in Figure 4(a)). These two values were affected by anomalies because in these two years the groundwater values were corrupted due to extraordinary pumping caused from a severe drought that had affected agriculture in 1988–1989 and in 2001–2002 (Hanel et al., 2018; Massarutto, 2000). A highly significant linear relationship, $y = a + b \times (F$-test $p \sim 0.00$, $R^2 = 0.81$), was obtained between the actual (y) and predicted (x) data, with the following parameter values in Equation (2): $A = 10.275$ (scaling parameter, converting the term in brackets to mm yr$^{-1}$), $B = 162.720$ cm (baseline groundwater depth, representing the 10th percentile) and $\eta = 0.5$ (empirical parameter). With the intercept $a = -0.053$ cm (±12 cm standard error) and the slope $b = 1.000$ (±0.084 standard error), only negligible departures of the data points from the 1:1 identity line are observed (Figure 4(a)).

310

315

320

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

**Figure 4.** Groundwater depth model evaluation. (a) Scatterplot of groundwater depth, Equation (2), in the calibration phase (period 1975–2012) at Cadriano. The inner bounds show confidence limits with a probability of 0.90 (hot pink area) and the outer bounds show prediction limits with a probability of 0.95 for new observations (light pink). The relative histogram of residuals is shown in the upper left corner, along with their normal fit (blue line). Two orange data points are outliers that were discarded from the calibration (see text). (b) Extended evaluation phase (period 1977–1988) comparing modelled (orange line) and observed (black line) groundwater depth (cm) for piezometer T12-Tonale (Giussano, Milan) during the period 1977–1998. (c) Similar to (b) referring to groundwater depth (m) for piezometer T16-Tortona during the period 2002–2017.
<figure><img src="image_16.png" /><figcaption>**Table 1.** Summary statistics of actual and modelled data (calibration).</figcaption></figure>

<table><thead><tr><th>Statistics</th><th>Actual data (calibration data set)</th><th>Model estimates, Equation (2)</th></tr></thead><tbody><tr><td>Count</td><td>38</td><td>38</td></tr><tr><td>Mean (cm)</td><td>138.1</td><td>138.1</td></tr><tr><td>Standard deviation (cm)</td><td>23.1</td><td>20.7</td></tr><tr><td>Coefficient of variation (%)</td><td>-17</td><td>-15</td></tr><tr><td>Minimum (cm)</td><td>-183.8</td><td>-177.6</td></tr><tr><td>Maximum (cm)</td><td>-89.3</td><td>-103.3</td></tr><tr><td>Range (cm)</td><td>-94.5</td><td>-74.3</td></tr><tr><td>Standardized skewness</td><td>-0.031</td><td>-0.051</td></tr><tr><td>Standardized kurtosis</td><td>-0.792</td><td>-1.170</td></tr></tbody></table>

According to the Durbin-Watson statistic (DW = 1.78, *p* = 0.22), there is no evidence of serial autocorrelation in the residuals of the calibration time series. The mean absolute error, MAE = 8.1 cm, is smaller than the standard error of the estimates (10.3 cm), and the mean absolute percentage error, MAPE = 6.1%, is below the performance target of excellence (MAPE < 10%). The corner histogram in Figure 4(a) shows a skew-free (Gaussian-like) distribution of the model residuals. As the maximum distance between the cumulative distributions of the two samples (actual and predicted) is equal to Dn = 0.11, the two-sided large sample Kolmogorov-Smirnov test (KS = 0.47, *p* = 0.98) suggests that the two data samples most likely come from the same distribution.

The calibration statistics with standardized skewness and kurtosis of the observed and estimated time series are within the expected Gaussian distribution range of -2 and +2 (Table 1). The highest *p*-value of the independent input variables is ~0.00, belonging to North Atlantic Oscillation + Atlantic Multidecadal Variability. Consequently, the GREM cannot be simplified.

## Extended model evaluation

An extended evaluation was performed using experimental measurements of the water table depth from the piezometric station T12-Tonale (Giussano, Milan), for the period 1977-1998, and the piezometric station T16-Tortona, for the period 2002-2017, which emphasized the need for parameter recalibration (e.g., adjusting parameters A and B of Equation (2) for variations in groundwater depth and the size of the area near the well). In fact, a shallower water table (as Cadriano) is mainly and rapidly affected by rainfall in the local area, whereas a deeper water table (as T12-Tonale) maintains the capacity to collect water over a longer time and from a wider variety of sources in the basin (e.g., from the mountain ridge). These hypotheses are empirically reflected in the calibrated parameter values: A = 20.5 and B = -240 cm at Tortona, A = 500 and B = -8000 at Tonale. The recalibrated model accurately reproduced the observed groundwater depth values at both sites (Figure 4(b,c)). The correlation between actual and modelled data was statistically significant for both sites (F-test, *p* ~ 0.00), with R² = 0.55 for Tortona and R² = 0.73 for Tonale. While there was evidence of serial autocorrelation in the residuals for T12-Tonale (DW = 0.83; *p* ~ 0.00), there was none for T16-Tortona (DW = 2.47; *p* = 0.79). However, no significant differences were found between the actual and modelled distributions for either station (T12-Tonale: Dn = 0.18, KS = 0.60, *p* = 0.86; T16-Tortona: Dn = 0.19, KS = 0.53, *p* = 0.94).
<figure><img src="image_17.png" /></figure>

We acknowledge the difficulty in separating natural and anthropogenic influences on groundwater levels, especially in heavily abstracted drinking water aquifers. However, our model successfully reproduces historical trends in both shallow and deep aquifers, even far from the calibration site. While further testing is needed, our results support the operational use of the GREM to estimate past groundwater dynamics in the Eastern Po Plain.

360

## Historical evolution of hydrological drought and groundwater depth

Groundwater levels in the Po Plain are influenced by both rainfall and flooding, as described in historical records. The Italian engineer Giuseppe Cadolini (1805–1858) highlighted this in *Il progetto del Canale Emiliano*, showing that flood events significantly replenish groundwater resources (Cadolini, 1893, p. 727, Supplementary file). Contemporary studies support this, indicating that both rainfall and floods shape groundwater dynamics in various landscapes (Cervi & Tazioli, 2021; Mancini et al., 2020; Pilli et al., 2012). On the contrary, periods with limited flooding or hydrological drought align with lower groundwater levels, as shown in Figure 5.

365

370

Pink bands in Figure 5 reflect drought periods identified through departures from the 11-year moving window groundwater Gaussian filter mean (bold grey horizontal line) of –147.3 cm (±19.4 cm standard error, Figure 5(c)). These bands indicate frequent droughts during the MCA and the MWP phases compared to the LIA. This further highlights that floods and groundwater levels are closely linked. Low groundwater levels can lead to fewer flood events and, eventually, droughts. To accurately describe ground-water drought, we need to set a specific threshold. For example, a drought occurs when groundwater levels fall below a certain level, such as the 25th percentile (in our case, equal to –164 cm), and ends when it rises again (Hanel et al., 2018). Such periods have occurred since the end of the MCA, with notable droughts in the years 1054–1105 and 1144–1165. In 1158–1159, Lancellotti (1680) recorded a 10-month drought (from May to the end of March) with a groundwater depth of –170 cm. Another major groundwater drought occurred between the years 1237 and 1253, with the 25th percentile reaching –174 cm. For 1243–1244, Lichtenthal (1854) described exceptionally dry conditions, with groundwater depth around –175 cm (about 29 cm below mean level; see also Figure 5(d)). The *Chronicon Parmense* (1038–1338) mentions the drought of 1237–1253, which was associated with the famine of 1247 (Bonazzi, 1902, pp. 672–673, Supplementary file).

375

380

385

The severe water shortage, which led to the famine, is also cited in climate reconstructions. Schoolman et al. (2018) used tree rings to show a drought in Italy between 1237 and 1253. Between 1100 and 1250, the MCA overlapped with the Great Solar Maximum, and its temperature peak occurred between 1000 and 1300, when a major temperature maximum occurred over most of Europe (Polovodova Asteman et al., 2018). Four mega-droughts occurred during the LIA: 1402–1447, 1610–1627, 1714–1740 and 1765–1806 (Figure 5, pink bands), which appear to be associated with below-average global temperatures (Evans et al., 2020). The first and last droughts roughly coincide with the megasequence of the ~1400–1480 and ~1770–1840 periods in central Europe (Ionita et al., 2021). The mid 15th-century mega-drought largely coincides with the Spörer Minimum (1420–1550), a period of low solar activity and major volcanic eruptions (Bauch, 2017). Famines in northern Italy between 1450 and 1559 likely resulted from climate trends and repeated crop failures

390

395

400
<figure><img src="image_18.png" /></figure>

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

**Figure 5.** Reconstructed evolution of groundwater depth in the period 1054–2020 CE. (a) Trasimeno lake anomaly level from 1401 to 2000 CE (Source: Burzigotti et al., 2003). (b) Low-extreme (25th percentile) groundwater depth (solid red line) with the corresponding linear trend (dashed red line). (c) Eleven-year Gaussian filter of the groundwater depth (black line), with mean (bold horizontal grey line) and standard deviation (grey bands), the latter calculated over an 11-year moving window (with each year centred on the fifth year of the window), and the main climatic periods with hydrological droughts (pink vertical bands). (d–g) Palmer Drought Severity Index for the years 1243–1244, 1615–1616, 1943–1945 and 2001–2020, respectively (data arranged via Climate Explorer, https://climexp.knmi.nl, based on the Old World Drought Atlas by Cook et al., 2015).
<figure><img src="image_20.png" /></figure>

(Alfani, 2010). Two long mega-droughts, 1721–1738 and 1765–1806, followed. The latter roughly corresponds to the Dalton Minimum (1790–1830), characterized by low solar activity, multiple volcanic eruptions and cooler temperatures (Ionita et al., 2021). Corradi (1972) describes this period as having dry, hot summers with several drought years. Ferrari (1723, p. 288, Supplementary file) specifically mentions the year 1718 as a severe one (groundwater depth equal to –170 cm).

405

Finzi (1986) also found four years of severe food shortages in Bologna: 1728, 1752, 1764 and 1766, the latter being the most disastrous after three years of drought. Other years with great hydrological drought were noted by Lichtenthal (1854) and Configliachi and Brugnatelli (1820): 1292 (groundwater depth equal to –165 mm), when the leaves of the trees dried up; 1334 (groundwater depth equal to –202 cm), when fields and vineyards burned due to heat and drought; 1391 (groundwater depth equal to –168 cm) and 1400–1401 (groundwater depth around –180 cm), when small rivers dried up. In 1403 (groundwater depth equal to –202 cm) and 1746 (groundwater depth equal to –171 cm), wells, rivers and vegetation dried up, as well as in 1540, which was so dry that green crops could not be planted (Targioni-Tozzetti, 1777). The extreme droughts of 1615 and 1616 (groundwater depth equal to –181 and –193 cm, respectively, as illustrated in Figure 5(e)), led to famine, followed by severe droughts in 1661 (Prosperi, 1989), 1665 (Anonymous, 1915), 1686 (Many authors, 1879), and 1775 and 1780, which were very hot and dry years (groundwater depth equal to –172 and –171 cm, respectively). One of the worst droughts on record occurred in 1763–1764. The severe water shortage, indicated by a groundwater depth of around –180 cm, led to a widespread famine in Italy in 1764 (Coppi, 1841). In addition, the year 1788 (groundwater depth equal to –193 cm), was marked by both extreme cold and persistent drought (SIASP, 1889, p. 445, Supplementary file).

410

415

420

425

The last climatic period, the MWP, evolved with successive drought waves: 1850–1893, 1930–1971 and 2006–2020. Major droughts affected France in 1893 and the Mediterranean in 1878 (Quintana Seguí et al., 2016). The late 19th century saw varying rainfall and drought in the Emilian Plain (Giusberti et al., 1999). The late 1930s and 1940s were particularly dry in Europe (Briffa et al., 1994). Northern Europe and Russia experienced severe droughts in the 1950s and 1960s, while Central Europe and the British Isles experienced droughts in the 1970s. The Mediterranean area and Baltic regions were affected in the 1990s and 2000s (Spinoni et al., 2015). Unlike previous mega-droughts, 21st-century droughts, such as 2003, 2015 and 2018 (Ionita et al., 2021), were aggravated by extreme heat. However, the 25th percentile of groundwater depth has decreased since 1606 (Figure 5(b)), with a decline of –0.85 cm per century (Mann–Kendall statistics S = 8471, Z = 3.11, *p* ~ 0.00), indicating a long-term drying trend in central Europe (Büntgen et al., 2021). The declining water level of Lake Trasimeno (43.13 N, 12.10 E, near Perugia, central Italy) since 1606 (Figure 5(a)) may offer insights into the Eastern Po Plain. However, caution is advised in directly extrapolating this downward trend to understand groundwater dynamics in this area, given the composition of the medium- to low-permeability rocks in the lake basin, coupled with extensive anthropogenic interventions over centuries, such as hydraulic works and catchment modifications (Burzigotti et al., 2003). The period 1965–2017 has seen increasing drought frequency and duration in the Po Plain, especially since 2001 (Figure 5(f)), together with declining river flows in the Mediterranean (García-Ruiz et al., 2011). Changes in precipitation patterns are a major driver of these trends (Baronetti et al., 2020).

430

435

440

445
# Conclusion

## Hydrological trends

This study provides a thorough historical analysis of groundwater depth trends in the Eastern Po Plain, synthesizing observational data, documentary evidence and proxy records to establish a long-term view of hydrological drought patterns. Our findings reveal that recent droughts, although severe, are not without historical precedent, highlighting a complex interplay between climatic variability and groundwater depletion. Observed declines in groundwater depth over time, particularly in the lower-percentile distributions, suggest that climate factors exert a stronger influence on groundwater levels that localized agricultural practices. While the site-specific measurements and estimates in this study cannot be extrapolated to the whole region with reasonable accuracy, they provide important insights into broader regional trends and highlight the role of climate as a primary driver of groundwater depletion.

Given the strong influence of climate variability on groundwater depletion, as evidenced by historical trends, there is a pressing need for adaptive and historically informed water management strategies in the Po Plain. By understanding how past climate patterns have shaped groundwater dynamics, policymakers can better anticipate and mitigate future risks

## Adaptive management strategies

This study highlights the need for adaptive and historically informed water management policies in the Po Plain. By linking historical climate variability to trends in groundwater depth, it provides a perspective on how past climate variability and associated groundwater dynamics have influenced water availability and, in extreme cases, caused socioeconomic stress. For example, we observe that historical droughts have often coincided with significant famines and agricultural challenges, such as those in the 13th and 18th centuries. Understanding these patterns allows water managers and policymakers to anticipate similar risks and implement more responsive water management strategies. This includes setting thresholds for groundwater extraction based on historical groundwater lows observed during mega-droughts, such as the 25th percentile groundwater depth of –174 cm reached during the major drought from 1237 to 1253, helping to prevent water scarcity escalating into socioeconomic crises during future droughts.

A key finding of the study is the important role that both flood events and droughts play in shaping groundwater dynamics. The dual influences of flood-induced recharge and drought-induced depletion suggest that policies should consider these cyclical groundwater behaviours. For instance, policies could allow temporary increases in groundwater extraction limits following significant flood events, leveraging the natural recharge of aquifers. Adaptive management that adjusts abstraction limits according to seasonal and long-term climate patterns would align water use with natural recharge processes, thereby reducing the risk of groundwater depletion during periods of reduced rainfall or limited flooding. In addition, by developing policy frameworks that integrate both aquifer-specific responses to climate events and broader climate drivers such as the North Atlantic Oscillation and the Atlantic Multidecadal Variability, regional water
<figure><img src="image_21.png" /></figure>

managers can better anticipate groundwater depletion risks during known climate phases.

490

To effectively address the increasing frequency and intensity of droughts, a shift from reactive to proactive water management is critical. Proactive steps could include setting lower abstraction quotas during drought years and planning for the use of alternative water sources, such as surface reservoirs, to relieve pressure on groundwater resources. In addition, the use of modern climate modelling techniques to simulate potential future scenarios could help to refine abstraction thresholds to prevent over-exploitation of groundwater during prolonged droughts, as observed historically during the Medieval Climate Anomaly and the Little Ice Age.

495

To effectively implement proactive water management strategies, it is essential to understand the large-scale climatic drivers that influence groundwater recharge. This study focuses on the North Atlantic Oscillation and Atlantic Multidecadal Variability as key indices that capture the atmospheric processes shaping European hydrology.

500

## Integration of large-scale climatic drivers

In assessing the robustness of the study, we acknowledge that alternative teleconnection indices, such as the Western Mediterranean Oscillation (WeMO) and the East Atlantic Pattern (EAP), could offer additional insights into regional climate variability. The WeMO, for instance, has been shown to significantly influence precipitation patterns in the western Mediterranean, particularly through its modulation of cyclonic activity and moisture transport (López-Moreno & Vicente-Serrano, 2008; Martín-Vide & López-Bustins, 2006). Similarly, the EAP, which describes variability in the North Atlantic jet stream, has been linked to shifts in atmospheric circulation patterns that affect temperature and precipitation across Europe (Comas-Bru & McDermott, 2014; Moore et al., 2013). However, the WeMO and EAP time series are only available from 1821 and 1950, respectively (via https://climexp.knmi.nl). Our decision to focus on the North Atlantic Oscillation and Atlantic Multidecadal Variability was motivated by their well-documented roles in modulating large-scale atmospheric processes that directly influence European hydrology, including precipitation patterns and temperature regimes (Hurrell, 1995; Sutton & Dong, 2012; Wang et al., 2017). These indices have been extensively validated and are widely used in hydrological studies, providing a robust framework for elucidating the climatic drivers of groundwater recharge in the Po Plain.

505

510

515

520

The North Atlantic Oscillation index used in this study was derived from instrumental sea-level pressure differences between the Azores High and the Icelandic Low, as described in Hurrell (1995). This approach is widely recognized for its high reliability in capturing the North Atlantic Oscillation's influence on European climate. For the Atlantic Multidecadal Variability index, we relied on sea-surface temperature anomalies in the North Atlantic, as reconstructed by Wang et al. (2017). These anomalies were derived from a combination of instrumental records and proxy-based reconstructions, ensuring consistency over the study period. While the North Atlantic Oscillation index is based on highly reliable instrumental SLP records, its use is restricted to periods with robust data to minimize uncertainties. The Atlantic Multidecadal Variability index, incorporating proxy-based reconstructions, introduces some uncertainty, particularly in earlier periods. However, the reconstruction by Wang et al. (2017) has been extensively

525

530
validated and is widely regarded as one of the most reliable Atlantic Multidecadal Variability data sets available.

We acknowledge that potential errors in the selected North Atlantic Oscillation and Atlantic Multidecadal Variability values – due to differences in calculation methods (e.g., station-based versus principal component approaches, as discussed in Pokorná & Huth, 2015) or differences in proxy-based reconstructions – could propagate into the groundwater depth estimates. To assess the robustness of our results, we tested alternative calculation methods, including variations in reference periods, spatial domains and methodological approaches. While the detailed results of these tests are not presented, preliminary findings confirmed the consistency of our results. In addition, the consistency of our findings across multiple assessment metrics, such as high goodness-of-fit (R² values ranging from 0.55 to 0.81) and distributional agreement supported by Kolmogorov-Smirnov tests, provides compelling evidence of the robustness and reliability of our chosen methodology. Future work may benefit from a systematic sensitivity analysis to further quantify the impact of alternative North Atlantic Oscillation and Atlantic Multidecadal Variability formulations on groundwater depth estimates, thereby strengthening the overall reliability of the results.

Although WeMO is known to affect precipitation patterns in the Mediterranean region (Català et al., 2005; Polcher et al., 2006), its regional focus may not fully capture the broader climatic forces relevant to our analysis. Similarly, while the EAP, which characterizes variability in the North Atlantic jet stream, could yield complementary insights, its influence overlaps considerably with that of the North Atlantic Oscillation in shaping European climate. Moreover, the North Atlantic Oscillation and Atlantic Multidecadal Variability indices used in this study were derived from standardized methods with well-documented uncertainties (Cook et al., 2015; Wang et al., 2017). Future research may benefit from exploring additional indices, such as WeMO or EAP, to further dissect local-scale climatic effects. However, for the purposes of this work, the North Atlantic Oscillation and Atlantic Multidecadal Variability indices provide the most comprehensive representation of the large-scale drivers influencing groundwater recharge in southern Europe. Integrating these insights into actionable water management policies will require continued collaborative efforts among researchers, policy-makers and local stakeholders.

## Research and policy integration

The clear link between climate cycles and groundwater droughts calls for integrated research that combines historical data with ongoing monitoring to ensure a dynamic approach to water resource management. Cross-sectoral collaboration between policymakers, climatologists and hydrologists can help to further refine groundwater depth thresholds and abstraction limits. Extending this research framework to other aquifer systems in the Po Plain and similar regions would not only strengthen the robustness of these thresholds, but also support climate-resilient, sustainable water management in different landscape types. While this study focuses on the Eastern Po Plain, extending this research framework to other aquifer systems across similar landscapes would not only strengthen the robustness of these thresholds but also support climate-resilient, sustainable water management strategies tailored to different environmental conditions.
<figure><img src="image_22.png" /></figure>

By integrating historical knowledge into contemporary water management practices, this research promotes a long-term view that balances immediate water needs with future availability. This approach could serve as a model for regions facing similar climate challenges, ultimately supporting a resilient and balanced water management strategy for the Eastern Po Plain and beyond.

In conclusion, this study highlights the critical role of historical climate variability and large-scale atmospheric forcing in shaping groundwater dynamics in the Eastern Po Plain. By integrating these findings into adaptive management strategies and fostering cross-sectoral collaboration, we can develop more resilient and sustainable approaches to water resource management in the face of increasing climate uncertainty.

## Acknowledgements

This research was performed as an investigator-driven study, without financial support. The data generated and analysed in this study are available in a data file published along with this article (Table S1.xls). The authors thank Prof. Francesca Ventura of DISTAL – Department of Agricultural and Food Sciences (University of Bologna, Italy) and the Experimental Farm of the University of Bologna (AUB) for providing the piezometric data of the Cadriano station used in this study.

## Disclosure statement

Q4 No potential conflict of interest was reported by the author(s).

## Authors contribution

N. D. developed the original research design and analysed the data. N. D., F. T. and G. B. wrote the article together and made the interpretations together. All authors reviewed the final manuscript.

## Q5 References

Alfani, G. (2010). Climate, population and famine in Northern Italy: General tendencies and Malthusian crisis, ca. 1450–1800. *Annales de démographie historique, 120*, 23–53. https://doi.org/10.3917/adh.120.0023

Anonymous. (1915). *Monumenti storici pubblicati dalla Regia Deputazione veneta di storia patria. Seria 3, Cronache e diarii, IV: Giovanni Soranzo*. Regia Deputazione Veneta di Storia Patria.

Baronetti, A., González-Hidalgo, J. C., Vicente-Serrano, S. M., Acquaotta, F., & Fratianni, S. (2020). A weekly spatio-temporal distribution of drought events over the Po Plain (North Italy) in the last five decades. *International Journal of Climatology, 40*(10), 4463–4476. https://doi.org/10.1002/joc.6467

Bauch, M. (2017). The day the sun turned blue: A volcanic eruption in the early 1460s and its possible climatic impact – A globally perceived volcanic disaster in the Late Middle Ages? In G. J. Schenk (Ed.), *Historical disaster experiences: A comparative and transcultural survey between Asia and Europe* (pp. 107–138). Springer. https://doi.org/10.1007/978-3-319-49163-9_6

Beven, K. J. (2012). *Rainfall-runoff modelling: The primer* (1st ed.). John Wiley & Sons. https://doi.org/10.1002/9781119951001

Bonazzi, M. (1902). *Chronicon Parmense ab anno 1038 usque ad annum 1338*. In L. A. Muratori (Ed.), *Rerum Italicarum Scriptores* (Vol. 9, pp. 672–673). Zanichelli.
<figure><img src="image_23.png" /></figure>

Briffa, K. R., Jones, P. D., & Hulme, M. (1994). Summer moisture variability across Europe, 1892–1991: An analysis based on the Palmer drought severity index. *International Journal of Climatology*, 14(5), 475–506. https://doi.org/10.1002/joc.3370140502

Büntgen, U., Urban, O., Krusic, P. J., Rybníček, M., Kolář, T., Kyncl, T., Ač, A., Koňasová, E., Čáslavský, J., Esper, J., Wagner, S., Saurer, M., Tegel, W., Dobrovolný, P., Cherubini, P., Reinig, F., & Trnka, M. (2021). Recent European drought extremes beyond Common Era background variability. *Nature Geoscience*, 14(4), 190–196. https://doi.org/10.1038/s41561-021-00698-0

Burzigotti, R., Dragoni, W., Evangelisti, C., & Gervasi, L. (2003, April 11–14). *The role of lake Trasimeno (central Italy) in the history of hydrology and water management* [Paper presentation]. IWHA 3rd International Conference, Alexandria, Egypt. https://www.researchgate.net/profile/Walter_Dragoni/publication/239587322_The_Role_of_Lake_Trasimeno_central_Italy_in_the_History_of_Hydrology_and_Water_Management/links/0deec52d99360da616000000/0deec52d99360da616000000.pdf

Casanueva, A., Rodríguez-Puebla, C., Frías, M. D., & González-Reviriego, N. (2014). Variability of extreme precipitation over Europe and its relationships with teleconnection patterns. *Hydrology and Earth System Sciences*, 18(2), 709–725. https://doi.org/10.5194/hess-18-709-2014

Català, M., Llasat, M. C., & Martín-Vide, J. (2005). The Western Mediterranean Oscillation and rainfall in the Iberian Peninsula. *International Journal of Climatology*, 25, 145–159. https://doi.org/10.1002/joc.1076

Cervi, F., & Tazioli, A. (2021). Quantifying streambed dispersion in an alluvial fan facing the northern Italian Apennines: Implications for groundwater management of vulnerable aquifers. *Hydrology*, 8(3), 118. https://doi.org/10.3390/hydrology8030118

Clifton, C., Evans, R., Hayes, S., Hirji, R., Puz, G., & Pizarro, C. (2010). *Water and climate change: Impacts on groundwater resources and adaptation options* (Water Working Notes No. 25). World Bank. http://hdl.handle.net/10986/27857

Comas-Bru, L., & McDermott, F. (2014). Impacts of the EA and SCA patterns on the European 20th-century NAO–winter climate relationship. *Quarterly Journal of the Royal Meteorological Society*, 140(679), 354–363. https://doi.org/10.1002/qj.2158

Configliachi, P., & Brugnatelli, G. (1820). *Giornale di fisica, chimica, storia naturale ed arti, Tomo III* [Journal of physics, chemistry, natural history, medicine and the arts, Volume 3]. Fusi.

Conte, M. (1992). La “bomba” meteorologica sul mediterraneo: Una climatologia sinottica [The meteorological “bomb” over the Mediterranean: A synoptic climatology]. *Bollettino Geofisico*, 25, 1–18.

Cook, E. R., Seager, R., Kushnir, Y., Briffa, K. R., Büntgen, U., Frank, D., Krusic, P. J., Tegel, W., van der Schrier, G., Andreu-Hayles, L., Baillie, M., Baittinger, C., Bleicher, N., Bonde, N., Brown, D., Carrer, M., Cooper, R., Čufar, K., Dittmar, C., ... Zang, C. (2015). Old World megadroughts and pluvials during the Common Era. *Science Advances*, 1(10), e1500561. https://doi.org/10.1126/sciadv.1500561

Coppi, A. (1841). *Discorso sull'agricoltura dell'Agro romano: Letto nell'Accademia Tiberina il dì 17 luglio 1837 da A. Coppi* [Speech on the agriculture of the Roman countryside read in the Tiber Academy on 17 July 1837]. Tipografia Salviucci.

Corradi, A. (1972). *Annali delle epidemie occorse in Italia dalle prime memorie fino al 1850* [Annals of epidemics that occurred in Italy from the first memories up to 1850]. Arnaldo Forni Editore.

Cramer, W., Guiot, J., Fader, M., Garrabou, J., Gattuso, J.-P., Iglesias, A., Lange, M. A., Lionello, P., Llasat, M. C., Paz, S., Peñuelas, J., Snoussi, M., Toreti, A., Tsimplis, M. N., & Xoplaki, E. (2018). Climate change and interconnected risks to sustainable development in the Mediterranean. *Nature Climate Change*, 8(11), 972–980. https://doi.org/10.1038/s41558-018-0299-2

Cuthbert, M. O., Taylor, R. G., Favreau, G., Todd, M. C., Shamsudduha, M., Villholth, K. G., MacDonald, A. M., Scanlon, B. R., Kotchoni, D. O. V., Vouillamoz, J.-M., Lawson, F. M. A., Adjomayi, P. A., Kashaigili, J., Seddon, D., Sorensen, J. P. R., Ebrahim, G. Y., Owor, M., Nyenje, P. M., Nazoumou, Y., & Kukuric, N. (2019). Observed controls on resilience of

620

625

630

635

640

645

650

655

660

665

670
<figure><img src="image_24.png" /></figure>

groundwater to climate variability in sub-Saharan Africa. *Nature*, 572(7768), 230–234. https://doi.org/10.1038/s41586-019-1441-7

De Lorenzo, G. (1939). Geologia e geografia fisica di Leonardo da Vinci [Geology and physical geography by Leonardo da Vinci]. In M. Dei Lavori Pubblici (Ed.), *Annali dei lavori pubblici [Annals of public works]* (pp. 233–263). Istituto Poligrafico dello Stato.

De Vita, P., Allocca, V., Manna, F., & Fabbrocino, S. (2012). Coupled decadal variability of the North Atlantic Oscillation, regional rainfall and karst spring discharges in the Campania region (southern Italy). *Hydrology and Earth System Sciences*, 21(5), 1389–1399. https://doi.org/10.5194/hess-16-1389-2012

Di Matteo, L., Spigarelli, A., & Ortenzi, S. (2021). Processes in the unsaturated zone by reliable soil water content estimation: Indications for soil water management from a sandy soil experimental field in Central Italy. *Sustainability*, 13(1), 227. https://doi.org/10.3390/su13010227

Diodato, N., Ljungqvist, F. C., & Bellocchi, G. (2020). Monthly storminess over the Po River Basin during the past millennium (800–2018 CE). *Environmental Research Communications*, 2(3), 031004. https://doi.org/10.1088/2515-7620/ab7ee9

Endo, T. (2022). The Kabu-ido system and factors affecting local groundwater extraction control: Case study of a customary groundwater management in Japan. *Water History*, 14(2), 169–191. https://doi.org/10.1007/s12685-022-00302-1

Esit, M., Kumar, S., Pandey, A., Lawrence, D. M., Rangwala, I., & Yeager, S. (2021). Seasonal to multi-year soil moisture drought forecasting. *Npj Climate and Atmospheric Science*, 4(1), 16. https://doi.org/10.1038/s41612-021-00172-z

Esposito, S., Beltrano, M. C., De Natale, F., Di Giuseppe, E., Iafrate, L., Libertà, A., Parisse, B., & Scaglione, M. (2015). *Atlante italiano del clima e dei cambiamenti climatici [Italian atlas of climate and climate change]*. Council for Agricultural Research and Agricultural Economy Analysis, Research unit for Climatology and Meteorology applied to Agriculture. https://www.reterurale.it/flex/cm/pages/ServeBLOB.php/L/IT/IDPagina/16319

Evans, S., Williams, G. P., Jones, N. L., Ames, D. P., & Nelson, E. J. (2020). Exploiting Earth observation data to impute groundwater level measurements with an extreme learning machine. *Remote Sensing*, 12(12), 2044. https://doi.org/10.3390/rs12122044

Fauzia, F., Surinaidu, L., Rahman, A., & Ahmed, S. (2021). Distributed groundwater recharge potentials assessment based on GIS model and its dynamics in the crystalline rocks of South India. *Scientific Reports*, 11(1), 11772. https://doi.org/10.1038/s41598-021-90898-w

Finzi, R. (1986). Il sole, la pioggia, il pane e il lavoro. Note su clima, raccolto, calendario agrario nel Bolognese durante il secolo xviii [The sun, the rain, the bread and the work. Notes on climate, harvest, agricultural calendar in the Bolognese area during the 17th century]. In R. Finzi (Ed.), *Le meteore e il frumento: Clima, agricoltura e meteorologia a Bologna nel 1900* [The meteors and wheat: climate, agriculture, meteorology in Bologna in the 18th century]. Il Mulino.

Garamhegyi, T., Kovács, J., Pongrácz, R., Tanos, P., & Hatvani, I. G. (2018). Investigation of the climate-driven periodicity of shallow groundwater level fluctuations in a Central-Eastern European agricultural region. *Hydrogeology Journal*, 26(3), 677–688. https://doi.org/10.1007/s10040-017-1665-2

García-Ruiz, J. M., López-Moreno, J. I., Vicente-Serrano, S. M., Lasanta-Martínez, T., & Beguería, S. (2011). Mediterranean water resources in a global change scenario. *Earth-Science Reviews*, 105, 121–139. https://doi.org/10.1016/j.earscirev.2011.01.006

Giusberti, F., Piro, F., & Sabbatani, S. (1999). *Acqua, ricchezza e salute: Il colera a Bologna nel XIX secolo* [Water, wealth and health: cholera in Bologna in the 19th century]. Editrice Compositori.

Hanel, M., Rakovec, O., Markonis, Y., Máca, P., Samaniego, L., Kyselý, J., & Kumar, R. (2018). Revisiting the recent European droughts from a long-term perspective. *Scientific Reports*, 8(1), 9499. https://doi.org/10.1038/s41598-018-27464-4

He, X., Bryant, B. P., Moran, T., Mach, K. J., Wei, Z., & Freyberg, D. L. (2021). Climate-informed hydrologic modeling and policy typology to guide managed aquifer recharge. *Science Advances*, 7(17), eabe6025. https://doi.org/10.1126/sciadv.abe6025

Hurrell, J. W. (1995). Decadal trends in the North Atlantic Oscillation: Regional temperatures and precipitation. *Science*, 269(5224), 676–679. https://doi.org/10.1126/science.269.5224.676

Q7

Q8
<figure><img src="image_25.png" /></figure>

Q9

Ionita, M., Dima, M., Nagavciuc, V., Scholz, P., & Lohmann, G. (2021). Past megadroughts in central Europe were longer, more severe and less warm than modern droughts. *Communications Earth & Environment, 2*(1), 61. https://doi.org/10.1038/s43247-021-00130-w

Kim, J. H., & Jackson, R. B. (2012). A global analysis of groundwater recharge for vegetation, climate, and soils. *Vadose Zone Journal, 11*(1), 188-222. https://doi.org/10.2136/vzj2011.0021RA

Kumar, C. P. (2012). Climate change and its impact on groundwater resources. *International Journal of Engineering Science, 1*, 43-60. http://www.researchinventy.com/v1i5.html

Kundzewicz, Z. W., & Döll, P. (2009). Will groundwater ease freshwater stress under climate change? *Hydrological Sciences Journal*, 54(4), 665-675. https://doi.org/10.1623/hysj.54.4.665

Lancellotti, D. S. (1680). *L'Hoggidi, overo il mondo non peggio ne più calamitoso del passato. Parte Prima [Today or the world is no worse nor more calamitous than the past. First Part]*. Gio Francesco Valuasense.

Lasagna, M., Mancini, S., De Luca, D. A., & Cravero, M. (2019). Piezometric levels in the Piedmont plain (NW Italy): Trend and hydrodynamic behaviour of the shallow aquifer. *Rendiconti Online della Società Geologica Italiana, 48*, 2-9. https://doi.org/10.3301/ROL.2019.30

Li, B., Rodell, M., Sheffield, J., Wood, E., & Sutanudjaja, E. (2019). Long-term, non-anthropogenic groundwater storage changes simulated by three global-scale hydrological models. *Scientific Reports*, 9(1), 10746. https://doi.org/10.1038/s41598-019-47219-z

Lichtenthal, P. (1854). *Manuale di geografia fisica. Seconda edizione originale rifusa e considerevolmente accresciuta* [Manual of physical geography. Second original edition recast and significantly enlarged]. Giuseppe Cioffi.

López-Moreno, J. I., & Vicente-Serrano, S. M. (2008). Positive and negative phases of the wintertime North Atlantic Oscillation and drought occurrence over Europe: A multitemporal-scale approach. *Journal of Climate, 21*(6), 1220-1243. https://doi.org/10.1175/2007JCLI1739.1

Mancini, C. P., Lollai, S., Volpi, E., & Fiori, A. (2020). Flood modeling and groundwater flooding in urbanized reclamation areas: The case of Rome (Italy). *Water, 12*(7), 2030. https://doi.org/10.3390/w12072030

Many authors. (1879). *Annali di Statistica, Serie 2, Volume 5* [Annals of Statistics, Series 2, Volume 5]. Tipografia Eredi Botti.

Martín-Vide, J., & López-Bustins, J. A. (2006). The Western Mediterranean Oscillation and rainfall in the Iberian Peninsula. *International Journal of Climatology, 26*(11), 1455-1475. https://doi.org/10.1002/joc.1388

Massarutto, A. (2000). *Agriculture, water resources and water policies in Italy*. University of Udine. https://doi.org/10.2139/ssrn.200151

McCallum, J. L., Dogramaci, S., Bai, A., Cook, P. G., Engdahl, N. B., Simmons, C. T., Skrzypek, G., & Grierson, P. F. (2020). Assessing temporal changes in groundwater recharge using spatial variations in groundwater ages. *Water Resources Research*, 56(8), e2020WR027240. https://doi.org/10.1029/2020WR027240

Meli, M., & Romagnoli, C. (2022). Evidence and implications of hydrological and climatic change in the Reno and Lamone river basins and related coastal areas (Emilia-Romagna, Northern Italy) over the last century. *Water, 14*(17), 2650. https://doi.org/10.3390/w14172650

Moore, G. W. K., Renfrew, I. A., & Pickart, R. S. (2013). Spatial and temporal variability of sea ice in the Labrador Sea and its impact on deep convection. *Journal of Physical Oceanography, 43*, 2099-2113. https://doi.org/10.1175/JPO-D-13-048.1

Neves, M. C., Jerez, S., & Trigo, R. M. (2019). The response of piezometric levels in Portugal to NAO, EA, and SCAND climate patterns. *Journal of Hydrology, 568*, 1105-1117. https://doi.org/10.1016/j.jhydrol.2018.11.054

Oke, I. A. (2007). Reliability and statistical assessment of methods for pipe network analysis. *Environmental Engineering and Sciences, 24*, 1481-1489. https://doi.org/10.1089/ees.2007.0001

Pilli, A., Sapigni, M., & Zuppi, G. M. (2012). Karstic and alluvial aquifers: A conceptual model for the plain - Prealps system (northeastern Italy). *Journal Hydrology, 464-465*, 94-106.

Q10

725

730

735

740

745

750

755

760

765

770

775
<figure><img src="image_26.png" /></figure>

Pokhrel, Y., Felfelani, F., Satoh, Y., Boulange, J., Burek, P., Gadeke, A., Gerten, D., Gosling, S. N., Grillakis, M., Gudmundsson, L., Hanasaki, N., Kim, H., Koutroulis, A., Liu, J., Papadimitriou, L., Schewe, J., Schmied, H. M., Stacke, T., Telteu, C.-E., ... Wada, Y. (2021). Global terrestrial water storage and drought severity under climate change. *Nature Climate Change*, 11(3), 226-233. https://doi.org/10.1038/s41558-020-00972-w

780

Pokorná, L., & Huth, R. (2015). Climate impacts of the NAO are sensitive to how the NAO is defined. *Theoretical and Applied Climatology*, 119(3-4), 639-652. https://doi.org/10.1007/s00704-014-1116-0

785

Polcher, J., Laval, K., & Picon, L. (2006). The impact of the Mediterranean Sea on the West African monsoon: A regional climate model study. *Quarterly Journal of the Royal Meteorological Society*, 132, 843-861. https://doi.org/10.1256/qj.05.23

Polovodova Asteman, I., Filipsson, H. L., & Nordberg, K. (2018). Tracing winter temperatures over the last two millennia using a north-east Atlantic coastal record. *Climate of the Past*, 14(7), 1097-1118. https://doi.org/10.5194/cp-14-1097-2018

790

Prosperi, A. (1989). *Storia di Cesena III. La dominazione pontificia* [History of Cesena III. Pontifical domination]. Bruno Ghigi Editore.

Quintana Seguí, P., Martin, E., Sánchez, E., Zribi, M., Vennetier, M., Vicente-Serrano, S. M., & Vidal, J. P. (2016). Drought: Observed trends, future projections. In S. Thiebault & J. P. Moatti (Eds.), *The Mediterranean under climate change* (pp. 123-131). IRD Editions. https://hal.archives-ouvertes.fr/hal-01401386

795

Rodell, M., Famiglietti, J. S., Wiese, D. N., Reager, J. T., Beaudoing, H. K., Landerer, F. W., & Lo, M. H. (2018). Emerging trends in global freshwater availability. *Nature*, 557(7707), 651-659. https://doi.org/10.1038/s41586-018-0123-1

800

Rust, W., Bloomfield, J., Cuthbert, M., Corstanje, R., & Holman, I. (2022). The importance of non-stationary multiannual periodicities in the North Atlantic Oscillation index for forecasting water resource drought. *Hydrology and Earth System Sciences*, 26(9), 2449-2467. https://doi.org/10.5194/hess-26-2449-2022

Saccò, M., Blyth, A. J., Humphreys, W. F., Cooper, S. J. B., White, N. E., Campbell, M., Mousavi-Derazmahalleh, M., Hua, Q., Mazumder, D., Smith, C., Griebler, C., & Grice, K. (2021). Rainfall as a trigger of ecological cascade effects in an Australian groundwater ecosystem. *Scientific Reports*, 11(1), 3694. https://doi.org/10.1038/s41598-021-83286-x

805

Schoolman, E., Mensing, S., & Piovesan, G. (2018). From the late medieval to early modern in the Rieti basin (AD 1325-1601): Paleoecological and historical approaches to a landscape in transition. *Historical Geography*, 46(1), 103-128. https://doi.org/10.1353/hgo.2018.0027

810

Sheng, H., & Chen, Y. Q. (2009, April 30-September 2). *Robustness analysis of the estimators for noisy long-range dependent time series* [Paper presentation]. ASME 2009 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, San Diego CA, USA. https://asmedigitalcollection.asme.org/IDETC-CIE/proceedings-abstract/IDETC-CIE2009/49019/1147/345706

815

Singh, N. K., & Borro, D. M. (2019). A Granger causality analysis of groundwater patterns over a half-century. *Scientific Reports*, 9(1), 12828. https://doi.org/10.1038/s41598-019-49278-8

Spinoni, J., Naumann, G., Vogt, J. V., & Barbosa, P. (2015). The biggest drought events in Europe from 1950 to 2012. *Journal of Hydrology: Regional Studies*, 3, 509-524. https://doi.org/10.1016/j.ejrh.2015.01.001

820

Štěpánek, P. (2005). *AnClim - software for time series analysis*. Institute of Geography, Faculty of Natural Sciences, Masaryk University.

Sutton, R. T., & Dong, B. (2012). Atlantic Ocean influence on a shift in European climate in the 1990s. *Nature Geoscience*, 5(11), 788-792. https://doi.org/10.1038/ngeo1595

825

Targioni-Tozzetti, G. (1777). *Relazioni d'alcuni viaggi fatti in diverse parti della Toscana: Per osservare le produzioni naturali, e gli antichi monumenti de essa, XI* [Reports of some trips made to different parts of Tuscany: To observe its natural productions and ancient monuments, XI]. Stamperia granducale.

Thornthwaite, C. W. (1948). An approach toward a rational classification of climate. *Geographical Review*, 38(1), 55-94. https://doi.org/10.2307/210739

830
Treppiedi, D., Cipolla, G., Francipane, A., & Noto, L. V. (2021). Detecting precipitation trend using a multiscale approach based on quantile regression over a Mediterranean area. *International Journal of Climatology*, 41(13), 5938-5955. https://doi.org/10.1002/joc.7161

Wang, J., Yang, B., Ljungqvist, F. C., Luterbacher, J., Osborn, T. J., Briffa, K. R., & Zorita, E. (2017). Internal and external forcing of multidecadal Atlantic climate variability over the past 1,200 years. *Nature Geoscience*, 10(7), 512-517. https://doi.org/10.1038/ngeo2962

Wu, W.-Y., Lo, M.-H., Wada, Y., Famiglietti, J. S., Reager, J. T., Yeh, P. J.-F., Ducharne, A., & Yang, Z.-L. (2020). Divergent effects of climate change on future groundwater availability in key mid-latitude aquifers. *Nature Communication*, 11(1), 3710. https://doi.org/10.1038/s41467-020-17581-y

Yang, Y., Guang, H., Batelaan, O., McVicar, T. R., Long, D., Piao, S., Liang, W., Liu, B., Jin, Z., & Simmons, C. T. (2016). Contrasting responses of water-use efficiency to drought across global terrestrial ecosystems. *Scientific Reports*, 6(6), 23284. https://doi.org/10.1038/srep23284

Zanchettin, D., Timmreck, C., Toohey, M., Jungclaus, J. H., Bittner, M., Lorenz, S. J., & Rubino, A. (2019). Clarifying the relative role of forcing uncertainties and initial condition unknowns in spreading the climate response to volcanic eruptions. *Geophysical Research Letters*, 46(3), 1602-1611. https://doi.org/10.1029/2018GL081018

Q14