# Assessing atmospheric influences for improving time-varying data-driven decadal predictions of Mediterranean spring discharge

Nazzareno Diodato, Francesco Granata, Fabio Di Nunno, Cristina Di Salvo, Gianni Bellocchi

► To cite this version:

Nazzareno Diodato, Francesco Granata, Fabio Di Nunno, Cristina Di Salvo, Gianni Bellocchi. Assessing atmospheric influences for improving time-varying data-driven decadal predictions of Mediterranean spring discharge. Hydrological Sciences Journal, 2025, pp.1 - 13. 10.1080/02626667.2025.2480128. hal-05098442

HAL Id: hal-05098442

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

Submitted on 5 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
<figure><img src="image_2.png" /></figure>

# Assessing atmospheric influences for improving time-varying data-driven decadal predictions of Mediterranean spring discharge

Nazzareno Diodato, Francesco Granata, Fabio Di Nunno, Cristina Di Salvo & Gianni Bellocchi

**To cite this article:** Nazzareno Diodato, Francesco Granata, Fabio Di Nunno, Cristina Di Salvo & Gianni Bellocchi (17 Mar 2025): Assessing atmospheric influences for improving time-varying data-driven decadal predictions of Mediterranean spring discharge, Hydrological Sciences Journal, DOI: 10.1080/02626667.2025.2480128

**To link to this article:** https://doi.org/10.1080/02626667.2025.2480128

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

View supplementary material

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

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

Accepted author version posted online: 17 Mar 2025.

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

Submit your article to this journal

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

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

View related articles

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

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

CrossMark

View Crossmark data

<figure><img src="image_11.png" /></figure>
**Publisher:** Taylor & Francis & IAHS

**Journal:** *Hydrological Sciences Journal*

**DOI:** 10.1080/02626667.2025.2480128

# Assessing atmospheric influences for improving time-varying data-driven decadal predictions of Mediterranean spring discharge

Nazzareno Diodato<sup>a</sup>, Francesco Granata<sup>b</sup>, Fabio Di Nunno<sup>b</sup>, Cristina Di

Salvo<sup>c</sup>, Gianni Bellocchi<sup>*a,d</sup>

<sup>a</sup> Met European Research Observatory – International Affiliates Program of the University Corporation for Atmospheric Research, Benevento, Italy.

ORCID: https://orcid.org/0000-0001-9549-1583

<sup>b</sup> Department of Civil and Mechanical Engineering, University of Cassino and Southern Lazio, Cassino, Italy.

ORCID: https://orcid.org/0000-0002-2268-6600

<sup>c</sup> Istituto di Geologia Ambientale e Geoingegneria, Consiglio Nazionale delle Ricerche, Montelibretti, Italy.

ORCID: https://orcid.org/0000-0002-8411-574X
d Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand,
France. ORCID: https://orcid.org/0000-0003-2712-7979

*Corresponding Author:* Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France. E-mail: **gianni.bellocchi@inrae.fr**
**Abstract**

Decadal climate forecasting for spring discharge is challenging due to its complex nature. Data-driven models (DDMs) are a promising approach. One such model, the time-varying Periodic AutoRegressive Moving Average with a Conditional Standard Deviation (PARMAX(TVAR)-CSD), applied to the Caposele spring (southern Italy) as a Mediterranean representative, exemplifies this approach by incorporating exogenous variables. Leveraging past data spanning 123 years (1899-2021) and encompassing both large-scale (Arctic Oscillation) and small-scale (precipitation) climate forcings, this model offers annual forecasts extending to 2060. Notably, a discernible interdecadal pattern emerged, projecting an increase in discharge during the decade 2031-2040, followed by a temporary drought phase from 2041 to 2050, consistent with global model projections. These findings highlight the influence of climate change and variability on precipitation patterns in southern Europe, with implications extending for global water resource management strategies, and emphasise the need for novel approaches to decadal forecasting of spring discharge in the Mediterranean.

**Keywords:** data-driven models; decadal climate forecasting; Mediterranean; spring discharge

**1. Introduction**

Groundwater, a critical resource for both human consumption and ecosystems (Fienen and Arshad, 2016), crucially contributing to achieving the Sustainable Development Goals (Figure 1), has faced increasing depletion across various regions worldwide, largely in response to climate variability (Alley et al., 2002; Dalin et al., 2017; Li et al., 2019; Rodell and Reager, 2023).

In particular, the decline in groundwater levels during the 20th century, driven by frequent and severe droughts (Chiang et al., 2021; Liu et al., 2023; Yang et al., 2023), has emerged as a major challenge for sustainable water management.
Accurately predicting spring discharge, which is highly sensitive to climate dynamics, is crucial for long-term resource planning, especially as demand rises (Raihan *et al.*, 2022). Despite its importance, long-term predictions of spring discharge remain challenging due to the complexity of the natural systems and the chaotic interplay between meteorological, hydrological and climate signals (Refsgaard and Storm, 1996). The difficulty is compounded by limited data and the increasing influence of anthropogenic activities on watersheds (Vu *et al.*, 2023). While deterministic models based on physical processes have been extensively used, their ability to capture the nonlinearity and uncertainty inherent in hydrological predictions, especially on decadal scales, remains limited (Bhattarai *et al.*, 2023). Global Circulation Models (GCMs) provide large-scale climate projections (Nerantzaki *et al.*, 2020; **Zhao et al., 2020**), but the uncertainties associated with random and systematic errors (Roy *et al.*, 2020; Huang and Zhao, 2022; Nolin *et al.*, 2023) often limit their accuracy in predicting local-scale hydrological processes, such as spring discharge. Deterministic forecasts of spring discharge are thus often unreliable, making data-driven models (DDMs) a promising alternative for hydrological forecasting.

In the context of groundwater prediction, DDMs - such as autoregressive
models - offer significant potential due to their ability to capture temporal patterns
in historical data (Aghelpour et al., 2022). Decadal-scale predictions are well
suited for implementation in autoregressive models, because they capitalise on
the predictability of long-term climate oscillations (Sun et al., 2018). These
models are particularly well-suited for predicting future values based on past
observations, with Autoregressive Integrated Moving Average (ARIMA) models
(Box and Jenkins, 1976) being widely adopted for time-series forecasting (Meehl et al., 2021; Vishwakarma et al., 2023).

DDMs approaches have been successfully used in forecasting sub-surface water flows, with studies using various approaches including ARIMA with external predictors (ARIMAX; Boggs et al. (2014), extreme value statistical decomposition (Zhong et al., 2016), ARIMA with seasonal component (SARIMA; Manakos et al., 2018), machine learning Hybrid Genetic Algorithm (Kovács and Stefanović, 2023), and long short-term memory Neural Network model (LSTM; Zhang et al. 2024). While these studies have considered various factors, they have neglected to include climate indexes as exogenous variables.

Studies on future hydroclimatic scenarios in Italy suggest a potential drought affecting spring discharge after 2030 (Diodato et al., 2017; Citrini et al., 2020). However, interdecadal predictions of groundwater flow remain exploratory due to limited knowledge and forecast skill for this timescale (Zhu et al., 2019). Addressing the limitations of time-varying factors in forecasting is ongoing. Given the sparsity of measurements for many hydrological variables (Gnann et al., 2023), representative sites becomes are crucial for coherent regional estimates. This study differs from previous work using non-linear autoregressive (NAR) models with exogenous input, which focussed on short-term predictions of spring flows in central Italy (Di Nunno et al., 2021). In recent decades, a decline in groundwater levels has been observed at the Caposele spring since 1980 (Fiorillo et al., 2021). In view of this observed trend, the present work takes on the challenge of predicting climate patterns for the coming decades.

To address this gap, our study introduces a novel Periodic AutoRegressive Moving Average with Conditional Standard Deviation - PARMAX(TVAR)-CSD -
model, specifically tailored to predict spring discharge on decadal scales. Established in
hydrology for capturing the periodic nature of time-series data (Vecchia, 1985; Bartolini
and Salas, 1993), Periodic Autoregressive Moving Average (PARMA) models (Ula,
1990) have seen further advancements with research on stochastic forecasting models
(Stagge and Moglen, 2013; Elganainy and Eldwer, 2018). Tesfaye et al. (2006)
identified periodic autoregressive moving average models for capturing seasonal
variations in river flows. Guo et al. (2023) introduced a time-varying stepwise
decomposition ensemble framework for nonstationary and nonlinear streamflow
forecasting on a monthly timescale. While these studies address the challenges of time-
varying factors, our work focuses specifically on decadal forecasting of spring
discharge, which is crucial for water resource management strategies. This model
uniquely incorporates exogenous mixed variables, with supporting input climate data
acting as soft information to mitigate forecast uncertainties. The inclusion of time-
varying (TVAR) parameters offers a distinct advantage in tracking time-varying
frequency, allowing the selection of the most reliable forecast amidst chaotic lapse
times (Lit, 2020). Exogenous factors effectively reduce nonlinearity and prevent shifts
to the autoregressive model in input-output data (Schär et al., 2024). Research has
shown that incorporating teleconnection patterns into predictive models yields notably
superior results compared to models that do not include such exogenous forcing,
extending the accuracy of forecasts several years into the future (Liu et al., 2023).

The primary objective of our modelling approach is to reduce uncertainty in
interannual variability by considering both large-scale factors, such as the Arctic
Oscillation (AO), for its established teleconnection patterns (Dünkeloh and Jacobeit,
2003), and smaller-scale climate forcings, in particular precipitation. The study also
recognises the dynamic shifts influenced by stochastic contributions from mid-latitude
atmospheric circulation and ocean-atmosphere connections, like the Atlantic Multidecadal Oscillation (Clement et al., 2015), which exerts a significant influence on atmospheric variability over the Mediterranean (Krichak et al., 2014).

## 2. Data and exploratory analysis

### 2.1 Study site

The geographical focus of this study is southern Italy (Figure 2a, b), characterised by a significant rainfall gradient (Figure 2b, c) owing to the convergence of the Apennine mountains and the southern Tyrrhenian coast. A key hydrological feature is the Caposele spring, which originates from the karst massif of Mount Cervialto (Figure 2d). The spring has an estimated recharge area of approximately 110 km² (Fiorillo et al., 2015), with a mean annual flow of about 4 m³ s⁻¹. The Caposele spring plays a crucial role in the regional water supply, particularly contributing to the *Acquedotto Pugliese* system in Apulia, southeastern Italy.

The region experiences a Mediterranean climate (Köppen-Geiger zone Cs), characterised by mild, wet winters and hot, dry summers (Beck et al., 2018). Variations in elevations result in climatic gradients, ranging higher temperatures at lower elevations Caposele (situated at 600 m a.s.l.) to cooler conditions at higher elevations (Cervialto Mountain at 1800 m a.s.l.).

### 2.2 Data sources and hydroclimatic characteristics

The CASD (Caposele Annual Spring Discharge) time-series data, spanning from 1899 to 2021 (Diodato et al., 2022), form the basis of this study. These data represent undisturbed conditions in the central Mediterranean, largely unaffected by human activity, making them valuable for hydrological modelling.
Exogenous input variables were derived from multiple sources: annual precipitation data for Montevergine (Diodato et al., 2023), and the Arctic Oscillation (AO) index from Trenberth and Paolino (1980) with updates available until 2021 from Climate Explorer (http://climexp.knmi.nl). The AO index is computed using Geopotential Height (GPH) data at 500 and 850 hPa levels, calculated as the zonal mean area-weighted 30° N-60° N mean minus the zonal mean area-weighted 60° N-90° N mean.

The integration of the AO index, indicative of large-scale climate variability, and the Montevergine precipitation, reflecting local precipitation, into an Atmospheric Pattern Index (section 3.3) provided insights into the hydroclimatic influences on the discharge dynamics of the Caposele spring.

## 2.3 Exploratory data analysis

To deeply explore and understand available data, we employed descriptive and exploratory data analysis involving visualising the autocorrelation function (ACF) and partial autocorrelation function (PACF), along with dynamic analyses in the space-time domain.

An initial analysis of the CASD data (Figure 3a) revealed a significant downward trend (Mann-Kendall test *p*<0.01). In addition, the variance of the CASD data was found to decrease over time (Mann-Kendall test *p*<0.05; Figure 3c). These findings indicate that the CASD time series is likely non-stationary in both first and second order.

Given the observation of varying variance and autocorrelation over time in the CASD data, we focused on models that could account for these characteristics. In addition, we considered heteroscedasticity in the residuals resulting from mean detrending (Figure 3b). The McLeod and Li test (1983) revealed no significant evidence
of heteroscedasticity in the residuals. This suggests that an autoregressive conditional heteroskedasticity (ARCH) model, typically used to handle heteroscedasticity, might not necessary for these specific data. Thus, we focused on finding a suitable $PARMAX(1,1)m \times CSD$ model to adequately fit the CASD data.

Figure 4a and Figure 4b illustrate the estimated ACF, showcasing the correlations between CASD levels at different time delays over the residual values. In this context, *lag k* represents the correlation between the standardised residuals of spring discharge at time *t* and time *t* + *k*, while taking into consideration all lower delays. This information is particularly valuable for determining the appropriate order of the autoregressive model required to best suit the data. Additionally, it is worth noting the presence of 95% confidence limits near the origin (0). If the probability bounds of that lag do not contain the calculated coefficient, there is a statistically significant association at the 95% confidence level. The ACF shows a linearly decreasing pattern, while the PACF decreases exponentially, with the first two spikes in the ACF occurring at lags 1 and 3.

Following Valenzuela *et al.* (2008), when the PACF drops more abruptly than the ACF, and the number of coefficients above the error criterion in the PACF is fewer than in the ACF, the model can be identified as ARMA(1,1), with an AR(*P*) component where *P*=1, and an MA(*Q*) component, where *Q* equals the number of PACF coefficients above the error criterion (in this case, 1) in its initial stage.

We turned to the space-time attractor to assess forecast predictability. An attractor is a set toward which a dynamically perturbed system converges after a sufficiently long period. The forecast can be deemed reliable if the system evolves in an organised manner after disturbance; otherwise, it is considered unreliable. The reliability parameters of the spring discharge time-series forecast were then determined, which is a critical step in DDMs. In the time-state domain, the CASD attractors, which represent the preferred
evolutionary path for the system, exhibited a higher degree of entanglement around a
directional branching (Figure 4c, red ellipse), indicating a potential predictability of
patterns, assuming accurate projections can be generated.

**3. Methods and model development**

**3.1 Model development and parameter specification**

To model the non-stationary CASD time series data, we employed a
PARMAX(TVAR)-CSD model, which integrates autoregressive components and
exogenous variables. The selection of the best-performing model involved careful
calibration of autoregressive parameters for the CASD time series. The
PARMAX(TVAR)-CSD model was chosen due to its effectiveness in detecting non-
stationary signals and using information from exogenous variables, with the final model
featuring an autoregressive component *P*=0 and a moving average component *Q*=2, based
on the validation results.

We employed a multi-metric approach to assess model performance. Root mean
square error (RMSE), mean absolute error (MAE) and mean absolute percentage error
(MAPE) were used to quantify error magnitudes. In addition, Pearson's correlation (r)
and Kling-Gupta efficiency (KGE) were used to assess the strength of the linear
relationship and overall agreement between predictions and observations. This multi
approach guided the selection of autoregressive parameters, ultimately maximising
model accuracy.

The Generalised Error Distribution (GED) was incorporated into the model to handle
non-stationary with fatter tails (McNeil *et al.*, 2015), accounting for extreme events like
intense rainfall events or drought (Cade and Noon, 2003). This distribution is well-
suited for forecasting spring discharge, where statistical properties like mean and
variance change over time (Vrugt et al., 2022). Additionally, the model’s autoregressive component includes *m* periodic elements, which must sum to zero. To meet this requirement, the first *m*−1 components are estimated, and their conjunction with the zero-sum constraint determines the final component (Lit et al., 2020).

## 3.2 Mathematical formulation

Figure 5 illustrates the conceptual framework presenting the essential phases of the PARMAX(TVAR)-CSD modelling approach used in this study.

In the first phase, we assume that the random variable *y* at time *t* (*y*<sub>*t*</sub>), depends not only on *y*<sub>*t*−*i*</sub>, but also varies with a constant period *m*. To represent the seasonal (periodical *m*) member in the seasonal autoregressive model, we employ an additive component, where the process is primarily made of two components.

The first is a *location* component (Montgomery *et al.*, 2008):

$$ \dot{y}_{t+1} = S_{t-i}(\text{ARX}) + m_{t-jm} \quad (1) $$

where *S*<sub>*t*</sub> represents the stochastic component, combining an autoregressive exogeneous (ARX) process with a deterministic segment *m*<sub>*t*</sub> of seasonal duration *m*. This formulation allows *y*<sub>*t*</sub> to be regarded as a process with predictable periodic behaviour (Montgomery *et al.*, 2008).

When observations are considered to be drawn from a probability distribution *p*(*y*<sub>*t*</sub> | *μ*<sub>*t*</sub>, *σ*<sub>*t*</sub>) with time-varying location *μ*<sub>*t*</sub> and scale *σ*<sub>*t*</sub> for *t* = 1,...,T (where T is the length of the timeseries), the location and scale are updated over time using the scaled score of the likelihood function. Generalised autoregressive score models, introduced by Creal *et al.* (2013) and adapted by Harvey (2013) for time-varying location/scale volatility models, provide a unified framework for incorporating time-varying
parameters into various nonlinear models (Lit et al., 2022). The link function $f(\cdot)$ relates the distribution parameters $\mu_t$ and $\sigma_t$:

$$ \mu_t = f(Z_t^\mu \dot{y}_t), \qquad \sigma_t = f(Z_t^\sigma \dot{y}_t) \tag{2} $$

where the selection matrices $Z_t^\mu$ and $Z_t^\sigma$ select components from the time-varying time scaled score vector $\dot{y}_t$. The time-varying parameter vector $\dot{y}_t$ is updated over time $\dot{y}_{t+1}$ and incorporates the stochastic deterministic components of a periodical scaled score time-varying autoregressive model (Lit, 2020, to which we refer for further insights):

$$ \dot{y}_{t+1} = \sum_{i=1}^{q} A_i \cdot s_{t-i+1} + \sum_{j=1}^{p} B_j \cdot \dot{y}_{t-jm+1} + \sum_{i=1}^{p} \sum_{j=1}^{q} (A_i \cdot B_j t_{t-i-j+1}) \cdot (\dot{y}_{t-i-jm}) + X \cdot \dot{y}_{t-i} + \omega_t \tag{3} $$

Here, $A$ and $B$ are fixed coefficient matrices, $s_t$ is the scaled score function driving the updating equation, $X$ represents the exogenous variable, and $\omega$ is an unknown coefficient representing a noise uncorrelated process with zero mean and constant variance $\sigma^2$. $A$ and $B$ depend on the static parameter vector $\psi$, which is estimated using Maximum Likelihood.

The definition of $s_t$ is:

$$ s_t = S_t \cdot V_t, \tag{4} $$

$$ \text{where: } V_t = \frac{\partial \log p(y_t | \dot{y}_t, \mathcal{F}_{t-1}; \psi)}{\partial \dot{y}_t} \tag{5} $$

for $t = 1,...,T$, and $V_t$ is the density score vector $p(y_t | \dot{y}_t, \mathcal{F}_{t-1})$. The scaling method considers the variance of the score and inverse Fisher information for score scaling. The information set $\mathcal{F}_{t-1}$ often includes lagged variables of $\dot{y}_t$ and $y_t$, as well as exogenous variables.
The second *log-scale* component is:

$$ \sigma = \exp(\text{Level} + X \cdot \beta + \text{Score}(1)) \tag{6} $$

The Score(1) is a similar factor as in Eq. (6), which lies in the ability of score-driven models to estimate the best time-varying parameter according to a Kullback-Leibler, i.e. an optimal filtering mechanism with respect to the true data generation *process* (Li *et al.*, 2022).

### 3.3 Integration of exogenous variables

In this study, we harnessed the power of the PARMAX(TVAR)-CSD to extract and assess the direct influence of exogenous climatic variables within a combined framework. The Atmospheric Pattern Index (API), a key indicator of climate variability in the central Mediterranean, served as our exogenous data. This index encompasses both small-scale (precipitation) and large-scale (Arctic Oscillation) climate influences. This intricate interplay between groundwater and climate change has been a topic of considerable interest among hydrogeologists (Klove *et al.*, 2014; Huo *et al.*, 2016). Consistent with these findings, earlier studies (Peleg *et al.*, 2012; Diodato and Fiorillo, 2013) highlighted the profound influence of the Mediterranean climate, characterised by substantial variability in annual and multiannual precipitation, on spring discharge. We introduced the API into the model components specifically designed for location and log-scale variations, with X in Eqs. (3) and (6), respectively.

To forecast CASD time-series, we incorporated nonlinear exogenous climate inputs through the API. Building on prior understanding of the natural temporal variaitions in CASD, unaffected by anthropogenic influence (Fiorillo and Guadagno, 2012), we
selected the PARMAX(TVAR)-CSD model. This model captures not only the time-series characteristics of natural variability in spring discharge but also the influence of exogenous climate forcing.

Our approach uses readily available annual precipitation (P) and annual mean Arctic Oscillation (AO) data. These features in the PARMAX(TVAR)-CSD model enable the capture of nonlinear behavior in an autoregressive time-series. In particular, the exogenous variable X in Eqs. (3) and (6) incorporates the API, which represents a combination of modified precipitation (MP) and Arctic Oscillation (AO) data. The API is defined as the ratio between MP at Montevergine station (located at an elevation of 1280 m a.s.l. and 40 km away from the Caposele aquifer) and AO, calculated as follows:

$$API = \frac{MP(MBOBS)}{(3+AO)} \qquad (7)$$

The modified precipitation, MP(MVOBS), is defined as:

$$MP(MVOBS) = \sqrt{P_{y=0}} + 0.6 \cdot \sqrt{P_{y=-1}} \qquad (8)$$

Eq. (8) considers the delayed effect of precipitation, representing the time it takes for water, after storage in the groundwater, to resurface through its hydrogeological cycle and contribute to spring discharge. The terms $P_{y=0}$ (current year) and $P_{y=-1}$ (previous year) in the equation reflect the precipitation amounts at these respective time points, capturing their influence on modified precipitation. The placement of AO at the denominator of the API formula reflects its negative correlation with spring discharge.

The API comprises two distinct components that jointly provide supportive information for forecasting CASD, operating at different spatial scales to characterise changes in spring discharge. Notably, in our analysis, when AO is considered independently, it exhibits an important correlation with variations in CASD output
(r=0.52). On the other hand, MP(*MVOBS*) exhibits a somewhat lower correlation with CASD (r=0.48). We assessed whether the integrated Eq. (7) provided a superior input for the PARMAX(TVAR)-CSD model compared to considering AO and PM(*MVOBS*) individually. The highest *p*-value among the independent variables was 0.000194, corresponding to MP(*MVOBS*). Since this *p*-value is less than 0.05, the term is statistically significant. Consequently, there was no need to remove any variables from the model. Thus, Eq. (7) serves as a suitable foundation to support the forecast of CASD.

The initial multiple regression model included the four independent variables MP(*MVOBS*), AO, NAO and AMO. To simplify the model, we iteratively removed variables with the highest p-values greater than 0.05. First, NAO was removed due to its p-value of 0.94. Then, AMO was excluded with a p-value of 0.65. The remaining model, with MR(*MVOBS*) and AO, showed all terms as statistically significant, with MR(*MVOBS*) having a p-value approaching 0. Thus, further simplification was unnecessary as all remaining terms were significant.

**3.4 Training, validation and forecasting phases**

The PARMAX(TVAR)-CSD model underwent a rigorous training phase using data
from 1899 to 1991 and was validated with data from 1992 to 2021. The validated model
was then applied to produce forecasts for the extended time horizon of 2022 to 2060. It
is important to note that there is no fixed ratio between the number of years for training
and validation, and the number of years of the forecasting horizon (Zhang *et al.*, 2022).
While general guidelines suggest that the training dataset should comprise more than
half of the whole dataset (Guo *et al.*, 2023), the specific characteristics of the spring
discharge time series in our study exhibit chaotic dynamics with notable predictability.
The use of the PARMAX(TVAR)-CSD model is considered appropriate for predicting spring discharge at Caposele spring, given the observed variability. Due to the unique dynamics of the dataset, a wider time window was used to compare the model predictions with climate simulation estimates. These climate simulations typically project into the mid-21st century (around 2050) and extend towards the end of the century (around 2100), e.g. Faggian et al. (2024). The CASD dataset, spanning the period from 1899 to 2021, owes its comprehensive and reliable nature to the efforts of Fiorillo and Guadagno (2010), who initiated the compilation of the time series, subsequently updated by Diodato et al. (2022).

Forecasts were generated using the Time Series Lab (Lit et al., 2022), Score Edition software V. 1.5 (https://timeserieslab.com), supported by Visual Recurrence Analysis (https://visual-recurrence-analysis.software.informer.com/4.9) and CurveExpert Professional 1.6 (https://www.curveexpert.net).

## 3.5 Assessment of decadal forecasts using CMIP5 models

We used decadal estimates from the Coupled Model Intercomparison Project Phase 5 (CMIP5, Du et al., 2022) to assess the consistency of our spring discharge forecasts (trends and interannual variability) with projected climate scenarios. Specifically, we focused on the RCP4.5 scenario, which represents a moderate greenhouse gas emission pathway (Moss et al., 2010). The Antecedent Annual Maximum of Total Precipitation on a Wet Day (AMTPWD) variable from CMIP5 models was selected because of its known relationship with spring discharge through delayed water return processes (Diodato et al., 2022). The visualisation of the forecasts involved plotting the predicted spring discharge time series alongside historical percentiles (15th and 95th) and the long-term mean discharge from the historical record.
## 3. Results and discussion

### 3.1 Performance assessment of PARMAX(TVAR)-CSD model

The visual representation in Figure 6a confirms a statistically significant association (p<0.05) between the API and yearly Caposele Annual Spring Discharge (CASD), supported by a strong Pearson's correlation coefficient (r=0.57, n=123, p<0.01). The majority of data points fall within the 95% prediction bounds, while only five data points lying outside of this limit. In addition, the Kolmogorov-Smirnov (K-S) test (p=0.52) reinforces the suitability of API as a factor supporting CASD forecasting, indicating no statistically significant difference between the distributions of CASD and API datasets.

This finding reaffirms the suitability of the API as a supportive factor for CASD forecasting. Hence, we incorporated the API index time series (2021-2060), combining modified precipitation and Arctic Oscillation data (provided in Appendix) into our spring discharge forecasting information stream.

Figure 6b depicts the close correspondence between observed (blue) and modelled (orange) data during the training (1899-1991) and the validation (1992-2021) periods using the PARMAX(TVAR)-CSD model. Error statistics for both the training and validation periods are presented in Table 1. The low error statistics (RMSE, MAE and MAPE) and high performance metrics (r and KGE) during validation demonstrate the model's effectiveness, and given the similarities between the validation (1992-2021) and the training (1899-1991) periods, it is reasonable to expect that the model will perform well in future projections. In addition, the model outperforms the PARMA(TVAR) model, which lacs exogenous variables, as shown by comparing the second and third columns of Table 1. Including an exogenous input improves the correct prediction of
spring discharge. The Log Loss function (0 to 1) measures uncertainty, but no set threshold exists. A MAPE of less than 10% in the PARMAX(TVAR)-CSD model indicates excellent accuracy (Osteragova and Osterag, 2012) for long-term CASD trend forecasts.

The robust performance of the PARMAX(TVAR)-CSD model across various metrics (RMSE, MAE, MAPE, r, KGE) demonstrates its accuracy and reliability. Crucially, the model exhibits low error (discrepancy between observed and estimated values), implying limited overall model uncertainty (Shrestha and Solomatine, 2008). This is likely due to the narrow parameter uncertainty (Lim *et al.*, 2006), as discussed in Diodato *et al*. (2012, 2018). While sensitivity analysis is a valuable tool for understanding how parameter variations affect model outputs, for a well-calibrated model with minimal error, such as PARMAX(TVAR)-CSD, such variations are likely to have minimal impact. This reduces the need for extensive sensitivity analysis in this case. This is consistent with the focus on model robustness in Sun *et al*. (2006) and the suggestion by Hyndman and Athanasopoulos (2020) to prioritise interpretation once robustness has been established. Similarly, Box and Jenkins (1976) acknowledge that sensitivity analysis may be redundant for consistently validated models with low error.

## 3.2 Spring discharge decadal forecast and implications for water resources planning

Understanding and forecasting spring discharge dynamics is a complex endeavour due to the interplay of temporal scales influencing groundwater processes. These scales encompass short-term variations and long-term trends, posing a significant challenge for managing spring groundwater systems effectively.

Figure 7 visualises annual spring discharge forecasts for the coming decades, indicating that CASD is expected to remain within the historical range from 1981 to
2020. These thick orange and blue lines represent the historical drought (15th percentile) and exceedance events (95th percentile), respectively, while the thin black line shows the long-term mean during this period. While the long-term trend does not reveal any discernible upward or downward trajectory, there are significant interdecadal oscillations, which are less pronounced than in the historical period. The forecasts suggest a general decline in spring discharge from the beginning of the forecast period, with notable interannual fluctuations around 2030.

Spring discharge levels are expected to hover around the long-term mean value during the current decade (2021-2030, Figure 7, thin black line). A notable shift is projected for the subsequent decade (2031-2040), with CASD surpassing the long-term mean value and potentially exceeding historical spring exceedance levels around 2036 (thick blue line). However, a sharp decrease in discharge is forecasted for the following decade (2041-2050). It is worth noting that only a few years are forecasted to approach the threshold value of spring drought (thick orange line).

The final decade (2051-2060) exhibits pronounced interannual oscillations, with transitions from drought-like conditions to periods of spring overshoot. These fluctuations suggest that changes in climate regimes could potentially jeopardize aquifer recharge and spring discharge rates.

A striking observation is the close alignment between the trend projected by our model and the Antecedent Annual Maximum of Total Precipitation on a Wet Day (AMTPWD) derived from 25 CMIP5 Models in the RCP4.5 scenario (Figure 7, dashed grey line). These CMIP5 models are: ACCESS1-0, bcc-csm1-1, bcc-csm1-1-m, CanESM2, CCSM4, CMCC-CM, CMCC-CMS, CNRM-CM5, CSIRO-Mk3-6-0, GFDL-ESM2G, GFDL-ESM2M, GISS-E2-R, GISS-E2-R, HadGEM2-CC, HadGEM2-ES, inmcm4, IPSL-CM5A-LR, IPSL-CM5B-LR, MIROC5, MIROC-ESM, MIROC-ESM-CHEM, MPI-ESM-LR,
MPI-ESM-MR, MRI-CGCM3 and NorESM1-M. AMTPWD, chosen for its stronger connection to spring discharge due to delayed water return (Diodato et al., 2022), exhibits a strong Pearson's correlation (r=0.60) with CASD. This consistency strengthens the link between groundwater flow and water intrusion (Costall et al., 2020) highlighting the importance of extreme rainfall events and AMTPWD as a proxy for spring discharge, as extreme rainfall events can lead to significant discharge and recharge cycles even in years with low precipitation. This phenomenon is observed not only in wet regions but also in dry areas worldwide (Cuthbert et al., 2019).

The projected decrease in spring discharge under the RCP4.5 scenario compared to the historical period is likely to be driven by two key factors. First, as documented by Diodato et al. (2023), precipitation in the mountainous areas surrounding the Caposele spring is projected to decrease significantly. Second, Diodato et al. (2025) indicate a concurrent and substantial reduction in future snowfall. Snowfall is critical for groundwater recharge, as it provides a slow and sustained release of water during the melting season that supports spring discharge. This expected decrease in snowfall, combined with the overall decrease in precipitation, suggests a reduced recharge capacity of the aquifer, ultimately leading to lower discharge rates in the future.

Our forecasting model shows consistent alignment with CMIP5 models, encompassing both the interannual and interdecadal oscillations. This alignment is particularly significant because it reflects the increasing climate variability observed worldwide, impacting spring discharge patterns (Kovács and Stevanović, 2023). This finding instills confidence in our forecasts and provides valuable insights for water resource planning in southern Italy, particularly where water resource management is challenging and regions are reliant on highly vulnerable aquifer systems to meet their water demands (Alfio et al., 2023). Accurate forecasts of CASD for upcoming decades
offer substantial advantages, particularly in operating the *Acquedotto Pugliese*. They facilitate strategic water supply management for Apulian communities, reducing shortages during low-flow periods. Furthermore, informed allocation and utilization of available water resources become possible. Accurate forecasts play are instrumental in enhancing preparedness for emergencies, particularly droughts. Early implementation of pre-emptive measures can minimise the potential impact on water supply. Furthermore, water infrastructure design can be adapted to accommodate long-term variations in spring discharge, contributing to overall system resilience. The knowledge gained from these forecasts empowers aqueduct managers to optimise supply system the operations, particularly pumping facilities. This optimisation not only ensures a steady water supply but also promotes sustainability and cost-effectiveness through reduced energy consumption.

The challenges posed by shifting climate patterns and their impact on spring discharge underscore the importance of proactive and adaptable water resource management strategies. While distinguishing between externally driven (natural and human induced) and internally generated decadal to multidecadal climate variability in complex models like GCMs remains challenging (Meehl et al., 2021), our approach provides valuable insights for time-scaling rainfall patterns that affect spring discharge.

### 3.3 Pitfalls and errors evaluating forecasts

Careful consideration must be given to potential pitfalls and errors when evaluating forecasts. Even well-developed models with low bias and variance can exhibit escalating prediction errors over time due to interdependencies between temporally close values (Haben, 2023). Non-linear autoregressive models with exogenous inputs,
despite their effectiveness, can be sensitive to initialisation effects. While these effects
may lead to an increase in prediction errors for successive years, the complexity is
compounded by variability across decades, with some periods exhibiting lower error
rates that potentially offset inaccuracies in long-term trends.

The reliance on input variables introduces another layer of sensitivity to errors. For
instance, discharge forecasting models heavily influenced by weather forecasts can be
inaccurate if the weather data itself are inaccurate due to measurement, forecasting or
calibration issues. Unexpected errors or inconsistencies in the data should prompt a
thorough examination of the input factors as potential sources of discrepancies.
External forcings acting on the model, such as the divergence between the rise in
precipitation and recharge relationship, introduce further sources of error. The
anticipated intensification of the hydrological cycle, resulting in more rainfall than
snowfall (Quante et al., 2021), adds to the complexity of forecast errors.
Acknowledging and addressing these complications is paramount to maintaining model
robustness and reliability. Future efforts should focus on refining data-driven models
while also comparing different model types.

**4. Conclusions**

**4.1 *Key findings and regional relevance***

This study highlights the effectiveness of a specific Autoregressive Moving Average
(ARMAX) model, PARMAX(TVAR)-CSD, in forecasting spring discharge in the
Mediterranean region. Our research has far-reaching implications for global water
resource management and climate adaptation, providing regional insights that transcend
geographical boundaries and offer valuable lessons. Accurate forecasting of spring
discharge enables proactive water management strategies in southern Italy. This
translates to optimising water infrastructure usage, ensuring reliable supplies during droughts and implementing timely water conservation measures. Ultimately, these forecasts inform infrastructure planning, enabling us to build more adaptable and resilient systems in the face of future changes in water availability.

With a changing climate and a growing population, effective management of water resources is critical. This research offers not only scientific advances, but also practical solutions to the global challenge of water security. The knowledge gained has transformative potential, particularly in managing the uncertainties associated with a changing climate. Our findings underline the need for further research into the complex relationship between climate variability and groundwater systems. This extends beyond southern Italy and contributes to broader climate adaptation strategies. The study promotes climate resilience by integrating climate models into local water management. This should enable policymakers to develop evidence-based policies aligned with global adaptation goals. It also paves the way for collaborative research on predictive modelling across regions.

## 4.2 Broader implications and future research directions

Sharing knowledge and best practices strengthens global adaptation efforts and leads to robust strategies for diverse climatic zones. Extending our methodology to other vulnerable aquifers can provide a more comprehensive understanding of the impacts of climate change on water resources. As access to clean water is a basic human right and essential for sustainable development, our research underlines the need for proactive and adaptive water management strategies.

The integration of advanced modelling techniques such as PARMAX(TVAR)-CSD model, coupled with knowledge of extreme weather events, provides a powerful tool to
address climate challenges. Improved understanding and forecasting capabilities will better equip us to ensure the availability of clean water for present and future generations. To strengthen the analysis and establish the relative merits of the PARMAX(TVAR)-CSD model, future research should include a benchmarking strategy that identifies the latest state-of-the-art scheme specifically designed for annual spring discharge forecasting and applies it to the Caposele spring discharge dataset, while also comparing it with other methods such as machine learning or specialised climate models, which would provide a broader context for understanding the strengths and limitations of PARMAX(TVAR)-CSD. Indeed, while recent studies have shown the effectiveness of deep learning approaches in karst spring discharge modeling (Wunsch et al., 2022), our study focusing on incorporating a time-varying autoregressive (TVAR) component allows the model to account for the changing relationships between past discharge and future values, potentially leading to more accurate forecasts (Ji and Ahn, 2023). Furthermore, this study relied on CMIP5 simulations for climate model assessment due to the unavailability of extreme variable outputs for CMIP6 (https://wcrp-cmip.org/cmip6) in the Climate Explorer platform at the time of analysis.

While CMIP6 represents an advancement in climate modelling, the daily maximum precipitation data required for our study had not yet been published. Future research should incorporate CMIP6 data as it becomes available, allowing for a more comprehensive evaluation of climate-driven hydrological changes and improving the robustness of predictive models. Moving forward, the inclusion of relevant climate indices, such as temperature or land use change, could be explored to further improve model accuracy. However, it is important to note that direct measurements of temperature and evapotranspiration are not currently available for the site in question. Even if these variables could be estimated indirectly by interpolation from global
datasets, their inclusion could introduce additional uncertainties in the prediction of spring discharge. Hence, any extension to include such variables would require careful selection based on their influence on spring discharge dynamics and their statistical compatibility with the existing framework. This approach promises to capture the complex dynamics of spring discharge under a changing climate with improved accuracy.

## Declarations

**Availability of data:** The set of raw data that support the findings of this study are available in the Supplementary material, Table S1.

**Funding:** The authors performed this research as an investigator-driven study with no financial support.

**Disclosure statement:** No potential conflict of interest was reported by the authors.

**Authors' contributions:** The contributions of the authors are as follows: N. D. developed the original research design, collected and analysed the historical documentary data. N. D., F. G., F. D. N., C.D.S. and G. B. contributed to the analysis and collaborated in writing the article and making interpretations. All authors participated in reviewing and finalising the manuscript.

## References

Aghelpour, P. et al., 2022. Comparing three types of data-driven models for monthly evapotranspiration prediction under heterogeneous climatic conditions. *Scientific Reports*, 12, 17363. https://doi.org/10.1038/s41598-022-22272-3

Alfio, M.R. et al., 2023. A comprehensive assessment of RCP4.5 projections and bias-correction techniques in a complex coastal karstic aquifer in the Mediterranean.
Frontiers in Earth Science, 11, 1231296.
https://doi.org/10.3389/feart.2023.1231296

Alley, W.M. et al., 2002. Flow and storage in groundwater systems. *Science*, 296, 1985–1990. https://doi.org/10.1126/science.1067123

Bartolini, P., Salas, J.D., 1993. Modeling of streamflow processes at different time scales. *Water Resources Research*, 29, 2573-2587.
https://doi.org/10.1029/93WR00747

Beck, H. E. et al., 2018. Present and future Köppen-Geiger climate classification maps at 1-km resolution. *Scientific Data*, 5, 1-12.
https://doi.org/10.1038/sdata.2018.214

Bhattarai, N. et al., 2023. Warming temperatures exacerbate groundwater depletion rates in India. *Science Advances*, 9, eadi1401 DOI: 10.1126/sciadv.adi1401

Blanc, E. et al., 2014. Modeling U.S. water resources under climate change. *Earth's Future*, 2, 197-224. https://doi.org/10.1002/2013EF000214

Boggs, K.G. et al., 2014. Forecast of natural aquifer discharge using a data-driven, statistical approach. *Ground Water*, 52, 853-863.
https://doi.org/10.1111/gwat.12133

Box, G.E.P. and Jenkins G., 1970. *Time series analysis: Forecasting and control*. San Francisco: Holden-Day.

Chiang, F., Mazdiyasni, O. and AghaKouchak, A., 2021. Evidence of anthropogenic impacts on global drought frequency, duration, and intensity. *Nature Communications*, 12, 2754. https://doi.org/10.1038/s41467-021-22314-w

Citrini, A., Camera, C. and Beretta, G.P., 2020. Nossana spring (Northern Italy) under climate change: Projections of future discharge rates and water availability. *Water*, 12, 387. https://doi.org/10.3390/w12020387
Clement, A., Bellomo K., Murphy LN, Cane MA, Mauritsen T., Rädel G., Stevens B., 2015. The Atlantic Multidecadal Oscillation without a role for ocean circulation. *Science*, **350**, 320-324. https://doi.org/10.1126/science.aab398

Costall, A.R. et al., 2020. Groundwater through flow and seawater intrusion in high quality coastal aquifers. *Scientific Reports*, **10**, 9866. https://doi.org/10.1038/s41598-020-66516-6

Creal, D.D., Koopman, S.J. and Lucas, A., 2013. Generalized autoregressive score models with applications. *Journal of Applied Econometrics*, **28**, 777-795.
https://doi.org/10.1002/jae.1279

Cuthbert, M.O. et al., 2019. Observed controls on resilience of groundwater to climate variability in sub-Saharan Africa. *Nature*, **572**, 230-234.
https://doi.org/10.1038/s41586-019-1441-7

Dalin, C. et al., 2017. Groundwater depletion embedded in international food trade. *Nature*, **543**, 700-704. https://doi.org/10.1038/nature21403

Diodato, N., Bellocchi, G., Bertolin, C., Camuffo, D., 2013. Mixed nonlinear regression for modelling historical temperatures in Central-Southern Italy. *Theoretical and Applied Climatology*, **113**, 187-196. https://doi.org/10.1007/s00704-012-0775-y

Diodato, N., Rubinetti, S., Bellocchi, G., 2025. Enhancing decadal snowfall forecasts in the Mediterranean mountains through informed atmospheric variability and climate data. *Hydrological Sciences Journal*, **70**, 294-310.
https://doi.org/10.1080/02626667.2024.2427356

Diodato, N., Fiorillo, F., 2013. Complexity-reduced in the hydroclimatological modelling of aquifer's discharge. *Water and Environment Journal*, **27**, 170-176.
http://onlinelibrary.wiley.com/doi/10.1111/j.1747-6593.2012.00336.x/abstract
Diodato, N. et al., 2017. Case study for investigating groundwater and the future of mountain spring discharges in Southern Italy. *Journal of Mountain Science*, 14, 1791-1800. https://link.springer.com/article/10.1007/s11629-017-4445-5

Diodato, N., Lanfredi, M., Bellocchi, G., 2023. Long-range, time-varying statistical prediction of annual precipitation in a Mediterranean remote site. *Environmental Research: Climate*, 2, 045011. https://doi.org/10.1088/2752-5295/acffe9

Diodato, N., Mao, L., Borrelli, P., Panagos, P., Fiorillo, F., Bellocchi, G., 2018. Climate-scale modelling of suspended sediment load in an Alpine catchment debris flow (Rio Cordon-northeastern Italy). *Geomorphology*, 309, 20-28. https://doi.org/10.1016/j geomorph.2018.02.026

Diodato, N. et al., 2022. Climatic fingerprint of spring discharge depletion in the southern Italian Apennines from 1601 to 2020 CE. *Environmental Research Communications*, 4, 125011. https://doi.org/10.1088/2515-7620/acae23

Di Nunno, F. et al., 2021. Prediction of spring flows using nonlinear autoregressive exogenous (NARX) neural network models. *Environmental Monitoring and Assessment*, 193, 350. https://doi.org/10.1007/s10661-021-09135-6

Du, Y., Wang, D., Zhu, J., Wang, D., Qi, X., Cai, J., 2022. Comprehensive assessment of CMIP5 and CMIP6 models in simulating and projecting precipitation over the global land. *International Journal of Climatology*, 42, 6859-6875. https://doi.org/10.1002/joc.7616

Dünseloh, A., Jacobeit, J., 2003. Circulation dynamics of Mediterranean precipitation variability 1948–98. *International Journal of Climatology*, 23, 1843-1866. https://doi.org/10.1002/joc.973
Elganainy, M.A., Eldwer, A.E., 2018. Stochastic forecasting models of the monthly streamflow for the Blue Nile at Eldiem Station. *Water Resources*, 45, 326-337.
https://doi.org/10.1134/S0097807818030041

Faggian, P., Trevisiol, A. and Decimi, G., 2024. Future projections of wet snow frequency and wet snow load on overhead high voltage conductors over Italy. *Cold Regions Science and Technology*, 217, 103980.
https://doi.org/10.1016/j.coldregions.2023.103980

Fienen M.N. and Arshad M., 2016. The International scale of the groundwater issue. In: A.J. Jakeman et al., ed. *Integrated groundwater management: Concepts, approaches and challenges*. Cham: Springer Nature, 21-48.

Fiorillo, F. and Guadagno, F.M., 2010. Karst spring discharges analysis in relation to drought periods, using the SPI. *Water Resources Management*, 24, 1867-1884.
https://doi.org/10.1007/s11269-009-9528-9

Fiorillo, F., 2011. The role of the evapotranspiration in the aquifer recharge processes of Mediterranean areas. In: G. Gerosa, ed. *Evapotranspiration – From measurements to agricultural and environmental applications*. Rijeka: InTech, 373-388.

Fiorillo, F. et al., 2021. Long-term trends in karst spring discharge and relation to climate factors and changes. *Hydrogeology Journal*, 29, 347-377.
https://doi.org/10.1007/s10040-020-02265-0

Fiorillo, F., Pagnozzi, M. and Ventafridda, G., 2015. A model to simulate recharge processes of karst massifs. *Hydrological Processes*, 29, 2301-2314.
https://doi.org/10.1002/hyp.10353
Gnann, S. et al., 2023. Functional relationships reveal differences in the water cycle representation of global water models. *Nature Water*, 1, 1079–1090. https://doi.org/10.1038/s44221-023-00160-y

Guo, B. et al., 2023. Prediction of groundwater level under the influence of groundwater exploitation using a data-driven method with the combination of time series analysis and long short-term memory: a case study of a coastal aquifer in Rizhao City, Northern China. *Frontiers in Environmental Science*, 11, 1253949. https://doi.org/10.3389/fenvs.2023.1253949

Haben, S., 2023. Verification and evaluation of load forecast models. In: S. Haben, M. Voss, W. Holderbaum, eds. *Core concepts and methods in load forecasting*. Cham: Springer International Publishing, 89-105. https://doi.org/10.1007/978-3-031-27852-5_7

Harvey, A.C., 2013. *Dynamic models for volatility and heavy tails: With applications to financial and economic time series*. Volume 52. Cambridge: Cambridge University Press.

Huang, Z. and Zhao, T., 2022. Predictive performance of ensemble hydroclimatic forecasts: Verification metrics, diagnostic plots and forecast attributes. *Wiley Interdisciplinary Reviews: Water*, 9, e1580. https://doi.org/10.1002/wat2.1580

Huo, A. et al., 2016. Groundwater storage and depletion trends in the Loess areas of China. *Environmental Earth Sciences*, 75, 1-11. https://doi.org/10.1007/s12665-016-5951-4

Hyndman, R.J., Athanasopoulos, G., 2020. Forecasting: principles and practice, 3rd edition. OTexts: Melbourne, Australia. https://otexts.com/fpp2

Ji, S., Ahn, K.-H., 2023. Temperature change-informed future multisite streamflow generation to support water supply vulnerability assessments under climate
change. *Journal of Hydrology*, 624, 129928.
https://doi.org/10.1016/j.jhydrol.2023.129928

Klove, B. et al., 2014. Climate change impacts on groundwater and dependent ecosystems. *Journal of Hydrology*, 518, 250–266.
https://doi.org/10.1016/j.jhydrol.2013.06.037

Kovács, A., Stevanović, Z., 2023. A combined stochastic-analytical method for the assessment of climate change impact on spring discharge. *Water*, 15, 629.
https://doi.org/10.3390/w15040629

Krichak, S.O. et al., 2014. Teleconnection–extreme precipitation relationships over the Mediterranean region. *Theoretical and Applied Climatology*, 117, 679–692.
https://doi.org/10.1007/s00704-013-1036-4

Li, B. et al., 2019. Long-term, non-anthropogenic groundwater storage changes simulated by three global-scale hydrological models. *Scientific Reports*, 9, 1074.
https://doi.org/10.1038/s41598-019-47219-z

Lim, K.J., Engel, B.A., Zang, T., Muthukrishnan, S., Choi, J., Kim, K., 2006. Effects of calibration on L-THIA GIS runoff and pollutant estimation. *Journal of Environmental Management*, 78, 35-43.
https://doi.org/10.1016/j.jenvman.2005.03.014

Lit, R., 2020. *Forecasting the VIX in the midst of COVID-19*. Time Series Lab - Article Series, https://timeserieslab.com/articles/rlit_vix_v2.pdf.

Lit, R., Koopman, S.J. and Harvey, A.C., 2022. *Time Series Lab, Score Edition*, https://timeserieslab.com.

Liu, K. et al., 2023. Past and future adverse response of terrestrial water storages to increased vegetation growth in drylands. *npj Climate and Atmospheric Science*, 6, 113. https://doi.org/10.1038/s41612-023-00437-9
Manakos, A. et al., 2019. Enhanced characterization of the Krania–Elassona structure and functioning allogenic karst aquifer in central Greece. *Geosciences*, 9 (1), 15. https://doi.org/10.3390/geosciences9010015

McLeod, A.I. and Li, W.K., 1983. Diagnostic checking ARMA time series models using squared residual autocorrelations. *Journal of Time Series Analysis*, 4, 269–273. https://doi.org/10.1111/j.1467-9892.1983.tb00373.x

Meehl, G.A. et al., 2021. Initialized Earth System prediction from subseasonal to decadal timescales. *Nature Reviews Earth & Environment*, 2, 340–357. https://doi.org/10.1038/s43017-021-00155-x

Nerantzaki, S.D., Hristopulos, D.T. and Nikolaidis, N.P., 2020. Estimation of the uncertainty of hydrologic predictions in a karstic Mediterranean watershed. *Science of the Total Environment*, 717, 137131. https://doi.org/10.1016/j.scitotenv.2020.137131

Nolin A.F. et al., 2023. Observed and projected trends in spring flood discharges for the Upper Harricana River, eastern boreal Canada. *Journal of Hydrology: Regional Studies*, 48, 101462, https://doi.org/10.1016/j.ejrh.2023.101462

Osteragova, E. and Osterag O., 2012. Forecasting using simple exponential smoothing method. *Acta Electrotechnica et Informatica*, 12, 62–66. https://doi.org/10.2478/v10198-012-0034-2

Peleg, N., Morin, E., Gvirtzman, H. and Enzel, Y., 2012. Rainfall, spring discharge and past human occupancy in the Eastern Mediterranean. *Climatic Change*, 112, 769–789. https://doi.org/10.1007/s10584-011-0232-4

Quante, L. et al., 2021. Regions of intensification of extreme snowfall under future warming. *Scientific Reports*, 11: 16621. https://doi.org/10.1038/s41598-021-95979-4
Raihan A.T., Bauer, S. and Mukhopadhaya, S., 2022. An AHP based approach to forecast groundwater level at potential recharge zones of Uckermark District, Brandenburg, Germany. *Scientific Reports*, 12, 6365. https://doi.org/10.1038/s41598-022-10403-9

Refsgaard, J.C. and Storm, B., 1996. Construction, calibration and validation of hydrological models. In: M.B. Abbott and J.C. Refsgaard, eds. *Distributed hydrological modeling*. Dordrecht: Kluwer Acad, 41–54.

Rodell, M. and Reager, J.T., 2023. Water cycle science enabled by the GRACE and GRACE-FO satellite missions. *Nature Water*, 1, 47–59. https://doi.org/10.1038/s44221-022-00005-0

Roy, T. et al., 2020. Global evaluation of seasonal precipitation and temperature forecasts from NMME. *Journal of Hydrometeorology*, 21, 2473–2486. https://doi.org/10.1175/jhm-d-19-0095.1

Schär, S., Marelli, S., Sudret, B., 2014. Emulating the dynamics of complex systems using autoregressive models on manifolds (mNARX). *Mechanical systems and signal processing*, 208, 110956. https://doi.org/10.1016/j.ymssp.2023.110956

Shrestha, D.L., Solomatine, D.P., 2008. Data-driven approaches for estimating uncertainty in rainfall-runoff modelling. *International Journal of River Basin Management*, 6, 109-122. https://doi.org/10.1080/15715124.2008.9635341

Stagge, J.H., Moglen, G.E., 2013. A nonparametric stochastic method for generating daily climate-adjusted streamflows. *Water Resources Research*, 49, 6179-6193. https://doi.org/10.1002/wrcr.20448

Sun, Y., Solomon, S., Dai, A., Portmann, R.W., 2006. How often does it rain? *Journal of Climate*, 19, 916–934. https://doi.org/10.1175/JCLI3672.1
Sun, Y. et al., 2018. Feasibility of predicting Vietnam's autumn rainfall regime based on the tree-ring record and decadal variability. *Climate, 6*, 42.
https://doi.org/10.3390/cli6020042

Trenberth, K. and Paolino, D., 1980. The Northern Hemisphere sea-level pressure data set: trends, errors and discontinuities. *Monthly Weather Review*, 108, 855-872.
https://doi.org/10.1175/1520-0493(1980)108<0855:TNHSLP>2.0.CO;2

Ula, T.A., 1990. Periodic covariance stationarity of multivariate periodic autoregressive moving average processes. *Water Resources Research*, 26, 855-861.
https://doi.org/10.1029/WR026i005p00855

Valenzuela, O. et al., 2008. Hybridization of intelligent techniques and ARIMA models for time series prediction. *Fuzzy Sets and Systems*, 159, 821-845.
https://doi.org/10.1016/j.fss.2007.11.003

Vecchia, A.V., 1985. Periodic Autoregressive-moving average (PARMA) modeling with applications to water resources. *JAWRA of the American Water Resources Association*, 21, 721-730. https://doi.org/10.1111/j.1752-1688.1985.tb00167.x

Vishwakarma, D.K. et al., 2023. Forecasting of stage-discharge in a non-perennial river using machine learning with gamma test. *Heliyon*, 9, e16290.
https://doi.org/10.1016/j.heliyon.2023.e16290

Vu, M.T. et al., 2023. Long-run forecasting surface and groundwater dynamics from intermittent observation data: An evaluation for 50 years. *Science of The Total Environment*, 880, 163338, https://doi.org/10.1016/j.scitotenv.2023.163338

Yang, Y. et al., 2023. Evapotranspiration on a greening Earth. *Nature Reviews Earth & Environment*, 4, 626-641. https://doi.org/10.1038/s43017-023-00464-3
Zhang, X.Y. et al., 2022. Advance prediction of coastal groundwater levels with temporal convolutional and long short-term memory networks. *Hydrology and Earth System Sciences*, 27, 83–96. https://doi.org/10.5194/hess-27-83-2023

Zhao, P. et al., 2020. Which precipitation forecasts to use? Deterministic versus coarser-resolution ensemble NWP models. *Quarterly Journal of the Royal Meteorological Society*, 147, 900–913. https://doi.org/10.1002/qj.3952

Zhong, Y. et al., 2016. A statistical model for karst spring discharge estimation under extensive groundwater development and extreme climate change. *Hydrological Sciences Journal*, 61, 2011-2023.

https://doi.org/10.1080/02626667.2015.1080825

Zhu, E., Yuan, X. and Wood, A.W., 2019. Benchmark decadal forecast skill for terrestrial water storage estimated by an elasticity framework. *Nature Communications*, 10, 1237. https://doi.org/10.1038/s41467-019-09245-3
# Appendix (Teleconnection-Climate Pattern Index forecasting)

## 1A. Initial modelling rules

When it comes to time-series forecasting methods, extrapolating future values from past data may differ from the approach used in the above CSD method. To understand the underlying behaviour of the time-series, it is essential to begin with a preliminary visual analysis of the data. Typically, this involves identifying cycles, analysing statistical data and looking for potential patterns like trends, cycles and oscillations around a mean, or even erratic and unpredictable variations. Additionally, it is important to examine fluctuations in the mean and variance, and consider the data's heteroscedasticity to establish stationarity.

To start an exploratory data examination, we initiated the development of the API time-series spanning from 1899 to 2021 (Appendix, Figure A1). The yearly API's evolution is displayed in Figure A1a, with the Mann-Kendall test (indicated by the grey dashed linear trend in Figure A1a) considerably supporting a decreasing trend (*p*<0.01). Figure A1b displays the residuals of this trend. While the API variance does not exhibit any apparent trend, there are some cyclical oscillations present (grey dashed linear trend in Figure A1c).

In summary, only the first order of the API time series appears to be stationary. This guides our study toward models that can account for heteroscedasticity simultaneously in the ARCH model, and where autocorrelation does not remain constant over time. Subsequently, we will evaluate the presence of heteroscedasticity using a test by examining the residual sequence (Figure A1b). Based on the McLeod-Li test results (*p*=0.16), heteroscedasticity is not present, and we should consider the appropriate PARMA(1,1)*m* model for matching the API data. Notably, the McLeod-Li test results also show a *p*=0.16.
## 2A. Model selection

For the reconstruction of the API for the period from 2022 to 2060, we adopted the PARMA(TVAR) model. The results of the McLeod-Li test indicate that the API time series is not heteroscedastic, and the validation of the PARMA(TVAR) model confirms that it performed effectively in this application. Consequently, we opted for this simplified approach within the PARMA(TVAR) model, which includes [Level + AR(2)+AR(1)+Periodical(25)+Score(1)]. In this model, we assume that the stochastic process $y(API)_t$ fluctuates with a constant period $m$ and is independent of $y(API)_{t-1}$. The seasonal (periodic) component of a seasonal autoregressive model is represented by an additive component, where the process is believed to be mainly influenced by a time-variable location component, as a log-scale component is assumed to be constant.

This location component follows the equation (Montgomery et al., 2008):

$$y(API)_t = S_{t-i}(AR) + m_{t-jm} \quad (A1)$$

where $S_t$ represents the stochastic component, which can be expressed as an autoregressive (AR) process combined with a deterministic member $m_t$ with the seasonal period of length $m$. According to Montgomery et al. (2008), $y(API)_t$ may be viewed in this context as a process with above-average predictability in its periodic behavior. The constant term is expressed as:

$$\log(\sigma_t^2) = K \cdot (S) \quad (A2)$$

## 3A. Parameterisation and evaluation of the PARMA(TVAR) model

In the Appendix, Figure A2 illustrates the progress of the observed and modelled data throughout the PARMA[1,1]25 training period (1899-1991) and the validation period (1992-2021). The closely aligned timelines in the validation period closely resemble those in the
training period, demonstrating the consistency between the observation (blue) and model (orange) lines. Error statistics for both the training and validation periods are shown in Table A1.

<table><caption>**Table A1.** PARMA(TVAR) model performance for the training (1899-1991) and validation (1992-2021) periods for the Rainfall Index time-series.</caption><thead><tr><th rowspan="2">Statistics</th><th colspan="2">PARMA(TVAR) Model</th></tr><tr><th>Training period (1899-1991)</th><th>Validation period (1992-2021)</th></tr></thead><tbody><tr><td>RMSE</td><td>3.9</td><td>4.2</td></tr><tr><td>MAE</td><td>3.1</td><td>3.4</td></tr><tr><td>MAPE (%)</td><td>12.3</td><td>16.2</td></tr><tr><td>Log Loss</td><td>4.2</td><td>4.0</td></tr></tbody></table>

For constructing model validation error statistics, the last 30 data values at the end of the time-series were used to. Confirming prediction inaccuracies is crucial to ensure the accuracy of the forecast. The first three statistics are used to assess the magnitude of errors. It is expected that the model will perform as predicted in future projection since the outcomes of the validation period (1992-2021) are comparable to those of the training period (1899-1991), as shown in Table A1.

Root Mean Square Error (RMSE) exhibits a particularly close resemblance between the calibration and validation states throughout the training phase. RMSE and MAE remain low during the validation state, with values of 4.2 and 3.4, respectively. Metrics like MAPE and Log Loss may underestimate the performance of the sub-
model, especially during calibration when the index is less than 15, as their focus on outliers or other data properties.

## Tables

<table><caption>**Table 1.** Performance comparison of PARMAX(TVAR)-CSD and PARMA(TVAR) models for Caposele spring discharge time-series during the training (1899-1991) and validation (1992-2021) periods. RMSE: root mean square error (optimal, 0; worst: +∞); MAE: mean absolute error (optimal, 0; worst: +∞); MAPE: mean absolute percentage error (optimal, 0; worst: +∞); Log Loss: logarithmic loss (optimum, 0; worst: +∞); r: Pearson's correlation coefficient (optimum, +1; worst: -1); KGE: Kling-Gupta efficiency (optimum, 1; worst: -∞).</caption><thead><tr><th rowspan="3">Statistics</th><th colspan="2">PARMAX(TVAR)-CSD Model</th><th>PARMA(TVAR) Model</th></tr><tr><th>Training period</th><th>Validation period</th><th>Validation period</th></tr><tr><th>(1899-1991)</th><th>(1992-2021)</th><th>(1992-2021)</th></tr></thead><tbody><tr><td>RMSE (m<sup>3</sup> s<sup>-1</sup>)</td><td>0.333</td><td>0.398</td><td>0.577</td></tr><tr><td>MAE (m<sup>3</sup> s<sup>-1</sup>)</td><td>0.243</td><td>0.326</td><td>0.458</td></tr><tr><td>MAPE (%)</td><td>6.0</td><td>9.11</td><td>12.7</td></tr><tr><td>Log Loss</td><td>0.25</td><td>0.75</td><td>2.31</td></tr><tr><td>r</td><td>0.73</td><td>0.56</td><td>0.49</td></tr><tr><td>KGE</td><td>0.71</td><td>0.56</td><td>0.30</td></tr></tbody></table>
## Figure captions

Figure 1. Illustration of the global water cycle and the central role of climate dynamics in achieving Sustainable Development Goals goals at different spatial scales. Water emerges as a critical factor in addressing the challenge of climate change and promoting socioeconomic growth. The narrative underscores the integral role of water as a key link between the climate system, the global economy and life-sustaining systems. Given recent significant changes in climate and increasing demand on water resources, advancing our understanding will require increasingly advanced approaches to predictive modelling (background image sourced from Freepik, http://www.freepik.com).

Figure 2. Environmental setting and hydrological framework of the Caposele-Cervialto aquifer. a. Geographical position within Europe, with Mediterranean outlined in the white square (arranged from https://www clipsafari com/clips/0188901-a-planetary-perspective-a-view-of-the-earth-from-space); b. Mean annual precipitation (1981-2010) for the Mediterranean and the broader study area (highlighted in red, arranged from: GPCC v2020 analysis (land) 0.25° resolution, via http://climexp knmi.nl); c. Downscaling of mean annual precipitation in the Campania region during 1971-2010, achieved through ordinary cokriging via ArcGIS Geostatistical Analyst - Environmental Systems Research Institute (data source: SCIA database, http://www.scia.isprambiente.it/wwwrootscia/Home_new.html); d. Perspective visualisation of the recharge area of the Caposele aquifer, with wooded shape in green, with reference to the Caposele station where discharge measurements are taken (landscape image from Google Earth and arranged with VanceAI OnlineTool).
Figure 3. Caposele spring discharge (CASD) data (1899-2021). a. Annual CASD time series data with regression line; b. Residuals of the annual CASD time series data resulting from mean detrending; c. Variance estimates of CASD time series data with an 11-year moving window and regression line.

Figure 4. Standardised Caposele spring discharge residuals and predictability (1899-2021). a. Autocorrelation function (ACF); b. Partial autocorrelation function (PACF), with horizontal black lines indicating the 95% confidence level boundaries; c. Time-series attractor in the phase-space domain (blue trajectory) displaying ellipse-like branching patterns (red contour).

Figure 5. Visual workflow for PARMAX(TVAR)-CSD model forecasting criteria (the picture in the corner at right is derived from; https://runoff.modelmywatershed.org).

Figure 6. Modelling process for forecasting Caposele spring discharge with an exogenous variable. a. Scatterplot illustrating the relationship between annual spring discharge from 1899 to 2021 and the exogenous API (Atmospheric Pattern Index). The inner (deep pink area) and outer (light pink area) bounds indicate the 90% and 95% confidence limits, respectively; b. Comparison of the modelled (orange line) and observed (blue line) annual spring discharge time series.

Figure 7. Time evolution of annual spring discharge, as observed from 1960 to 2021 (blue line) and as forecasted by Eq. (3) from 2022 to 2060 (red line), including related uncertainty bounds (grey band), long-term mean value (thin black horizontal line), and highlighting spring drought (thick orange horizontal line) and exceedance (thick blue line).
horizontal line). The dashed grey line represents the antecedent annual Maximum of
Total Precipitation on a Wet Day estimated as: *(MTPWD<sub>y</sub> + 0.75·MTPWD<sub>y-1</sub> + 0.50·MTPWD<sub>y-2</sub> - 0.75·MTPWD<sub>y-3</sub> - 0.50·MTPWD<sub>y-4</sub> - 0.25·MTPWD<sub>y-5</sub>)*, with
MTPWD forecasted by ensemble CMIP5 in RCP4.5 scenario, sourced from Climate
Explorer http://climexp.knmi.nl).

Figure A1. Atmospheric Pattern Index (API) data (1899-2021). a. Annual API time-series and its regression line; b. Its residuals; c. Variance estimates of API with an 11-year moving window and its regression line.

Figure A2 a. Modelled (orange line) and observed (blue line) annual API time-series during the validation stage; b. Modelled (orange line) and observed (blue line) annual Rainfall Index time-series during the forecasted stage.
<figure><img src="image_12.png" /><figcaption>Figure_1</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure_1A</figcaption></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>Figure-2</figcaption></figure>

<figure><img src="image_16.png" /></figure>
<figure><img src="image_17.png" /><figcaption>Figure_2A</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure_3</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure_4</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure_5</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure_6</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure-7</figcaption></figure>

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