# Neural Networks for Seasonal Groundwater Level Forecasting

Houa Bedjil, Guillaume Artigue, Jean Luchier, Rafik Ferraoun, Faustine

Bousquet, Anne Johannet

► To cite this version:

Houa Bedjil, Guillaume Artigue, Jean Luchier, Rafik Ferraoun, Faustine Bousquet, et al.. Neural Networks for Seasonal Groundwater Level Forecasting. Inżynieria Mineralna, 2025, 3 (2), 10.29227/IM-2025-02-03-24. [hal-05483340](https://hal-05483340)

HAL Id: hal-05483340

https://imt-mines-ales.hal.science/hal-05483340v1

Submitted on 29 Jan 2026

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 CC BY 4.0 - Attribution - International License
Abstract

Reliable and sustainable water resource management largely depends on the accurate prediction of groundwater levels, a challenge that has become increasingly critical due to climate change, and the growing exploitation of aquifers due to human activities. Reliable forecasts are essential for effective decision-making in water resource management, particularly in regions facing seasonal fluctuations or long-term declines in groundwater availability. The hidden nature of underground reservoirs makes it impossible to gain detailed knowledge of structural variability and water paths. Apart from large homogeneous basins, the physical properties of aquifers are thus poorly understood, making them difficult to model using physical or physically based models. Also, difficulties encountered in delimiting their feeding basin make it a challenge to fully understand them and accurately simulate and forecast their evolution. Thus, the lack of physical knowledge of underground hydrosystems encourages the use of systemic modeling methods, based on statistical approaches that do not require such knowledge. These include models based on artificial neural networks, particularly multilayer perceptrons (MLPs), that have proven to be a promising approach for improving the accuracy of hydrogeological forecasts. These models offer the advantage of capturing complex, nonlinear relationships between meteorological variables and groundwater level fluctuations, often outperforming traditional statistical methods. Although IA methods need a sufficient volume of data, they do not require extensive information about the hydrodynamics of the aquifer, making them easier to implement in poorly known contexts. In this study, we develop a seasonal groundwater level forecasting model for 30, 60, and 90 days ahead, for a piezometer located in the basement of the Sélune River watershed, in the Manche region of France. To achieve this, we use a multilayer perceptron (MLP) neural network while integrating regularization techniques to prevent overfitting and enhance generalization capabilities. The model is trained on a historical dataset covering a 15-year period with a decadal time step, incorporating several input variables, including temperature, actual evapotranspiration, precipitation, and the discharge of the nearby river. Cross validation and a separated test set are used to ensure the robustness of the model. Furthermore, the impact of future meteorological data on forecast quality is analyzed to assess its contribution to improving model performance. The results show a variation in forecasting accuracy depending on the lead time, with a general tendency toward reduced accuracy over longer periods. However, the integration of available future meteorological data significantly enhances model performance, highlighting its importance in hydrogeological forecasting and resource planning, opening new paths of inquiry for future developments. Estimating future trends in meteorological variables, even if uncertain, therefore appears to be a significant way of improving seasonal water level forecasts.

**Keywords:** water resource management, seasonal groundwater level forecasting, multilayer perceptron (MLP), future meteorological data

1. Introduction

A significant number of territories across Europe faced severe drought issues in the past few years. Climate change progressively started adding its effects to the natural variability, as well as the increasing consumption of water. In this context, water resources management becomes a key element of sustainability and questions our capacity to anticipate extreme events and to rely on a forecast of an accurate enough range of water tables. As the duration of droughts is several months or even years and the inertia of aquifers is often high, the lead time of such forecasts can be seasonal, i.e. from one to six months ahead.

Two main modelling approaches can be implemented for such purposes. First, it is possible to use physically based models, that attempt to represent the physical behaviour of the aquifer, from which it is theoretically possible to extract the water table at any spatial point. This approach requires i) a deep knowledge of the considered aquifer, which is time consuming, and ii) a forecast of the inputs to get forecasts of the outputs. In the context of an operational development on several aquifers, this is difficult to implement. For example, [1] combine GloSea5—a seasonal meteorological forecast model- with hydrogeological conceptual models and produce better forecasts than a persistence forecasting approach, although the contribution of GloSea5 to forecast quality is not always clearly identifiable.

<footnote><sup>1)</sup>HydroSciences Montpellier (Univ. Montpellier, IMT Mines Alès, IRD, CNRS), 6 Avenue de Clavières, 30319 Alès Cedex
France, houa-bedjil@imageau.com</footnote>

<footnote><sup>2)</sup> imageau - 158 Ave Dr Fleming 30900 Nimes France, houa-bedjil@imageau.com</footnote>

<footnote><sup>3)</sup> HydroSciences Montpellier (Univ. Montpellier, IMT Mines Alès, IRD, CNRS), 6 Avenue de Clavières, 30319 Alès Cedex
France, guillaume.artigue@mines-ales.fr, ORCID: 0000-0002-4610-584X</footnote>

<footnote><sup>4)</sup> imageau - 158 Ave Dr Fleming 30900 Nimes France, jean.luchier@imageau.com</footnote>

<footnote><sup>5)</sup> imageau - 158 Ave Dr Fleming 30900 Nimes France, rafik.ferraoun@imageau.com</footnote>

<footnote><sup>6)</sup> imageau - 158 Ave Dr Fleming 30900 Nimes France, faustine.bousquet@imageau.com</footnote>

<footnote><sup>7*)</sup> HydroSciences Montpellier (Univ. Montpellier, IMT Mines Alès, IRD, CNRS), 6 Avenue de Clavières, 30319 Alès Cedex
France, anne.johannet@mines-ales.fr; ORCID 0000-0001-5717-9871</footnote>
Second, it is possible to use statistical approaches, such as machine learning methods, that do not need extensive information about the physical behaviour and characteristics of the aquifer. They rely on a causal relationship between explicative inputs and desired outputs. Among these approaches, neural networks have proven their effectiveness producing forecasts of groundwater level without any forecast of their inputs (rainfall, temperature, evapotranspiration, etc.) on a short lead-time. In general, the lead time of these forecasts are limited to the response time of the system, which can be high for aquifers. For example, as early as 2005, [2] produces neural networks groundwater level forecasts of good apparent quality in Greece. However, the overexploitation of the aquifer, the good scores obtained are probably more due to the ability of the model to capture a general trend than to reproduce seasonal variability. More recent versions of such approaches are being developed: [3] used NARX models to forecast groundwater levels with a performance depending on the inertia of the considered systems.

Increasing the lead time of forecasts can be achieved by using meteorological seasonal forecasts. Unfortunately, despite the strong interest for this field, a review of the literature shows that it remains difficult to obtain accurate enough forecasts over one to two months, whether purely for forecasting meteorological variables [4, 5, 6, 7], or for precipitation forecasts aimed at extending the forecast horizon for groundwater levels [8, 9, 10]. The specificity stands in the fact that hybrid approaches exist, especially for downscaling or post-treatment of low resolution physically based models e.g [11, 12].

The north-western part of France is no exception to the above-described context. In this paper, we focus on the Manche area of France to provide 30, 60 and 90 days ahead forecasts of groundwater level, on an aquifer of interest for water resource management. The purpose is to allow the local managers optimizing water resource allocation, long time enough before a drought. As the selected approach should be used in operational conditions and relatively easy to generalize to other sites, we did not choose a physically based approach, but neural networks modelling. Besides, the availability of accurate enough meteorological seasonal forecasts is not guaranteed, enjoining us to assess the ability of the selected models to provide forecasts: (i) without any forecast of their inputs, or (ii) with very simple future input data. We thus implemented recurrent models for 30, 60 and 90 days ahead lead times with: (i) no future input data, (ii) perfect forecasts and (iii) forecasts of inputs based on past years' average. After presenting the site and the data, the model design and selection approach will be developed, before presenting and discussing the results obtained and the perspectives they open.

## 2. Materials and methods

### 2.1 Study Area: Basement of the Sélune Watershed

The Sélune watershed, covering an area of 1,014 km², extends across three French departments: Manche, Ille-et-Vilaine, and Mayenne. The Sélune originates at an altitude of 140 meters in Saint-Cyr-du-Bailleul, in the Manche department. With a length of 68 km, it flows into the Bay of Mont-Saint-Michel. Its average slope is 2‰ (Fig. 1).

Located on the Armorican Massif basement, the watershed is underlain by a geological structure composed of weathered and fractured rocks, resulting in mixed flow patterns: porous flow in weathered saprolites and predominantly fractured flow. Groundwater flows follow the topography, aligning hydrogeological and topographic watersheds.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Location and cartographic representation of the Sélune watershed, extent of the aquifer, measurement stations, elevation and hydrography.</figcaption></figure>

The data used comes from the Saint-Laurent-de-Terregatte piezometer, the Safran grid point provided by Météo-France (temperature, precipitation, and evapotranspiration) [13], and the hydrometric station of the Sélune in Saint-Aubin-de-Terregatte (table 1). This data covers the period from 2008 to 2022 with a decadal time step. It is worth noting that this basin exhibits strong interactions between groundwater and surface water [14].
Building on this dataset, the following figures provide a comprehensive overview of the hydroclimatic dynamics observed in the Sélune basin from 2008 to 2022. Fig. 2 displays the annual mean air temperature from 1991 to 2020, along with the 10 year's average for the periods 1991–2000, 2001–2010, and 2011–2020. The data show an overall warming trend, with the 10 year's average increasing from 11.53 °C in 1991–2000 to 12.05 °C in 2011–2020, illustrating the regional manifestation of climate change.

Throughout the year, precipitation is significant in all seasons, with higher amounts typically observed in autumn and winter due to the frequent oceanic fluxes from the Atlantic Ocean. Interannual variability is also strong, as some years -such as 2014- stand out as particularly wet, with totals exceeding 1,100 mm, while others, including 2010, 2016, and 2022, were notably drier, with cumulative precipitation occasionally falling below 800 mm.

<table><caption>Tab. 1. Overview of collected data.</caption><thead><tr><td>Variable</td><td>Data Source</td><td>Unit</td><td>Min value</td><td>Max value</td><td>Average</td></tr></thead><tbody><tr><td>Groundwater level</td><td>Saint-Laurent de Terregatte piezometer</td><td>m.a.s.l</td><td>44.84</td><td>47.48</td><td>45.65</td></tr><tr><td>Precipitation</td><td>Météo-France</td><td>mm</td><td>0</td><td>48.90</td><td>2.41</td></tr><tr><td>Temperature</td><td>Météo-France</td><td>°C</td><td>-4.3</td><td>31.10</td><td>11.38</td></tr><tr><td>Evapotranspiration</td><td>Météo-France</td><td>mm</td><td>-0.10</td><td>6.20</td><td>1.65</td></tr><tr><td>Discharge</td><td>Hydrometric site of the Sélune</td><td>m³/s</td><td>0.88</td><td>70.70</td><td>9.83</td></tr></tbody></table>

Fig. 2 also shows the evolution of river discharge and precipitation at 10-day time-step from 2008 to 2022. The hydrograph reveals strong seasonal and interannual variations, with discharge peaks generally occurring in winter and spring, in response to high precipitation. Extended low-flow periods during dry years highlight the sensitivity of streamflow to climatic conditions.

Fig. 3. presents the 10-day evolution of groundwater level, precipitation, and real evapotranspiration from 2008 to 2022. Groundwater levels exhibit strong seasonal dynamics, with peaks typically observed during late winter or early spring, corresponding to recharge periods driven by precipitation. Minimum levels occur during summer, in association with high evapotranspiration rates. Interannual variability is also evident in the amplitude of these seasonal fluctuations, reflecting the influence of varying climatic inputs. Years with limited winter precipitation or high summer evapotranspiration show reduced recharge and lower overall groundwater levels.

The joint analysis of climatic variables and hydrogeological responses reveals a close coupling between atmospheric forcing and water resources. Groundwater levels and streamflow are strongly dependent on the timing and quantity of precipitation, as well as on evapotranspiration rates.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of temperature (left), discharge and precipitation (right) during the 2008–2022 period.</figcaption></figure>

Years with combined low rainfall and high temperatures, such as 2017 and 2022, are especially critical, resulting in diminished recharge, increased aquifer depletion, and extended low-flow periods in rivers. These conditions intensify hydrological stress, threatening both ecosystems and water availability for human use.

These results emphasize the increasing vulnerability of water resources to climatic variability and underscore the importance of adaptive management strategies that integrate seasonal dynamics and long-term climate trends.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of groundwater level in response to precipitation and evapotranspiration (2008–2022).</figcaption></figure>
## 2.2 Neural networks models

A neuron is a mathematical operator that applies a generally non-linear transformation from its inputs to its outputs. In practice, neurons make a weighted sum of their inputs -called the neuron potential- and apply an activation function to it, to give an output.

Neural networks are composed of interconnected neurons organized into specific architectures. They consist of hidden neurons whose outputs are not the outputs of the model, and output neurons. One of these architectures is the multilayer perceptron (MLP) (Fig. 4), frequently used for its universal approximation properties, demonstrated by [15], stating that a MLP can approximate any differentiable function if it has enough hidden neurons and a high-quality dataset. [16] also demonstrates that MLP are more parsimonious than models linearly dependent on their parameters, as the number of parameters they use increases linearly with the number of variables. These properties make MLP models particularly interesting for modelling non-linear complex phenomena.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Multilayer Perceptron schematic representation.</figcaption></figure>

MLP can be classified based on how they manage time. Three types of models can be distinguished

* The static model: it only relies on exogenous inputs, with no dependence on previous measured or simulated outputs. Time has no functional role in these models (Eq. 1)

$$ \hat{y}(k) = \varphi(x(k), \dots, x(k - n_r + 1), \mathbf{W}) \qquad \text{Eq. 1} $$

Where $\hat{y}(k)$ is the estimated output at the discrete time $k$; $\varphi$ the nonlinear function implemented by the model; $\mathbf{x}$ the input vector; $n_r$ the sliding time windows size that defines the length of the necessary data; $\mathbf{W}$ the vector of parameters.

* The feedforward model: this model uses previous observed values of the modelled variable. One can assume that this additive data is used as a state variable for the model. The impact of this variable decreases as the lead time increases. (Eq. 2)

$$ \hat{y}(k) = \varphi(y(k-1), \dots, y(k-r); x(k), \dots, x(k-n_r+1), \mathbf{W}) \qquad \text{Eq. 2} $$

Where $y(k)$ is the observed value of the simulated variable at the discrete time $k$ and the other notations are the same than for the static model.

* The recurrent model: a model where the output depends on both exogenous inputs, and previous simulated outputs, creating a feedback loop (Eq. 3).

$$ \hat{y}(k) = \varphi(\hat{y}(k-1), \dots, \hat{y}(k-r); x(k), \dots, x(k-n_r+1), \mathbf{W}) \qquad \text{Eq. 3} $$

Where $r$ is the order of the recurrent model and the other notations are the same than for the static and the feedforward models.

According to [17], the use of feedforward models should be dedicated to the cases where the output noise is low and the state noise is high, and conversely for recurrent models. In our case, the output noise is low (high quality of outputs measurements), and the state noise is high (poor knowledge of the processes), but as the lead time is high, we assume that the impact of the previous observed outputs is low, especially for 60 and 90 days ahead. Recurrent models will thus be evaluated in this study.

## 2.3 Neural networks training, overtraining and regularization methods

Neural networks training consists in adjusting the parameters of the model to minimize a cost function. For this purpose, in this study, we use the Levenberg-Marquardt algorithm, with a least squares criterion as the cost function. The network parameters are initially assigned randomly and then optimized on a training set consisting of inputs-output vector pairs. The model adjusts its parameters based on performance evaluations until a stopping criterion is met. The primary goal is generalization, ensuring that the model performs well on unseen data, which is validated using a separate test set. A too complex model risks overfitting, reproducing the sum of the signal and of the associated noise, leading to a high variance of the modeled output and a low bias. Conversely, a too simple model may underfit, failing to capture important relationships and resulting in high bias and a low variance. This is called the “bias-variance dilemma” first described by [18]. To overcome this issue and achieve an optimal balance, regularization methods are used to constrain model complexity while maintaining its generalization capability.

First, early stopping [19] prevents overfitting by stopping the training process when the model starts to overtrains on the training data. This is achieved using a separate validation set, known as the stopping set. At each training iteration, the model's error is measured on the stopping set. When this error begins to increase, it indicates that the model is overtraining. Training is then halted, and the model parameters used are those of the model calculated at the stop iteration.

As the error calculated on the training set is not a relevant estimator of the generalization error [18], cross-validation [20] is applied to estimate the generalization error. The training set is thus divided into $V$ subsets, and the model is trained on $V-1$ subsets while being validated on the remaining one. This process is repeated $V$ times, ensuring that each subset serves as a validation set
once. The cross-validation score (S) is computed as the median of the validation scores across the V iterations. The model with the best cross-validation score is selected as the best-performing model.

As neural networks training is sensitive to the initial random values of the parameters, different initializations can lead to variations in the final model, affecting its performance and convergence. To mitigate this issue, several different initializations are performed to create an ensemble of models, offering a broader range of possible solutions and improving overall robustness [21, 22]. These three regularization methods are applied simultaneously.

## 2.4 Variable selection

The available data is divided into one-year-long subsets, spanning from January to December. These subsets are used for training (2008 to 2022), except for the ones allocated to the stop set (2021) and the test set (2016 and 2017). The stop set is selected as the most “average” year of the database, and the test set is selected as the most extreme years, here in terms of drought, which is the purpose of modelling.

Cross correlation analysis is used to select input variables by identifying those most correlated with the outputs. They are computed to determine the appropriate range of sliding input time windows investigation, thus limiting the number of combinations to evaluate during the model selection step. This analysis provides insights into the linear explanatory power of the inputs, the system's response time, its memory effect, and the reasonable lead time.

The analysis reveals the slow response dynamics of the aquifer system to input variations, highlighting its inertia. As a result, forecasts, without reliable meteorological prediction, remain theoretically feasible up to two months ahead, aligning with the system's response time (table 2).

<table><caption>Tab. 2. Maximum correlation, monthly estimates of response time (in time step of 10 days) and maximum correlation duration (in time step of 10 days) for each variable with respect to groundwater level.</caption><thead><tr><th>Variable</th><th>Cross-correlation indicator</th><th>Groundwater level</th></tr></thead><tbody><tr><td rowspan="3">Groundwater level</td><td>Maximum correlation value</td><td>1</td></tr><tr><td>Response time</td><td>-</td></tr><tr><td>Memory effect duration</td><td>2.5</td></tr><tr><td rowspan="3">Evapotranspiration</td><td>Maximum correlation value</td><td>0.71</td></tr><tr><td>Response time</td><td>4</td></tr><tr><td>Maximum duration of correlation</td><td>7</td></tr><tr><td rowspan="3">Temperature</td><td>Maximum correlation value</td><td>0.70</td></tr><tr><td>Response time</td><td>7</td></tr><tr><td>Maximum duration of correlation</td><td>9</td></tr><tr><td rowspan="3">Precipitation</td><td>Maximum correlation value</td><td>0.45</td></tr><tr><td>Response time</td><td>1</td></tr><tr><td>Maximum duration of correlation</td><td>4</td></tr><tr><td rowspan="3">Discharge</td><td>Maximum correlation value</td><td>0.93</td></tr><tr><td>Response time</td><td>0</td></tr><tr><td>Maximum duration of correlation</td><td>8</td></tr></tbody></table>

## 2.5 Quality criteria

To assess the reliability and performance of the forecasting models, two widely used statistical criteria are employed: the persistence criterion and the Nash-Sutcliffe efficiency (NSE). These metrics evaluate how well the predicted values reproduce the observed data, providing insight into the models' predictive skill and robustness.

*   Persistence Criterion

This criterion evaluates the deviation between the forecast and a naïve prediction. A prediction is considered naïve when the forecasted value at a lead-time $l$ ($k+l$) is simply equal to the observed value at time $k$, resulting in a persistence score of 0. A perfect forecast yields a score of 1 [23]. The persistence criterion is calculated as follows:

$$C_p = 1 - \frac{\sum_{k=1}^{n-l} (x(k+l) - y(k+l))^2}{\sum_{k=1}^{n-l} (x(k) - x(k+l))^2} \qquad \text{Eq. 4}$$

Where:

$l$ is the lead time, $k$ is the discrete current time, $x(k+l)$ is the observed value at the lead time, $y(k+l)$ is the predicted value at the lead time $l$, $n$ is the number of observed-forecast couples.

*   Nash-Sutcliffe Efficiency (NSE) or Coefficient of Determination (R²)

The Nash criterion [24] is widely used to assess the accuracy of forecasts. It quantifies how well the predicted values match observed values over the simulation period. Its value ranges from -∞ to 1. A score close to 1 indicates highly accurate predictions, while a value near 0 implies that the model performs no better than simply using the mean of the observed data. Negative values indicate that the model performs worse than the mean. According to [24], a model is considered satisfactory with a score of 0.7 and very good from 0.9 upwards. It is calculated as:

$$R^2 = 1 - \frac{\sum_{k=1}^{n-l} (x(k+l) - y(k+l))^2}{\sum_{k=1}^{n-l} (x(k) - \bar{x}(k))^2} \qquad \text{Eq. 5}$$

Where:

$l$ is the lead time, $k$ is the discrete current time, $x(k+l)$ is the observed value at the lead time, $y(k+l)$ is the predicted value at the lead time $l$, $n$ is the number of observed-forecast couples, $\bar{x}$ is the mean of the observed data.
## 2.6 Model design and selection

As previously stated, the model used in this study is recurrent. Both hidden neurons numbers and sliding time windows of inputs are explored using cross validation, to select the more robust models, for the three lead times and for the three kinds of meteorological variables forecasts used:

* Recurrent model (Rec), without any forecast of meteorological variables.

* Recurrent model, using perfect forecasts of meteorological variables (i.e. observed variables in the future) (Rec_perfect), which should be the best reachable performance,

* Recurrent model, using past averages of necessary variables as a forecast of these variables in the future (Rec_average).

By default, the model order is set to three, which is considered sufficient to capture the second derivative of the previously simulated outputs. The values of the average variables -evapotranspiration (ETR_avg), discharge (Q_avg), temperature (T_avg), and precipitation (RR_avg) are fixed up to the forecast horizon. Meanwhile, all other hyperparameters related to observed variables (denoted as 'obs') are explored, as detailed in Table 3. In the Rec_perfect configuration, it is assumed that the observed input variables ('obs') are available up to the forecast horizon (k + l), allowing the model to benefit from “ideal” future input conditions.

The general structure of the recurrent model, including both observed and average input variables, is illustrated in Fig. 5. This figure also highlights the recursive nature of the model, where three previous outputs (RecInputs) are used as part of the input for each prediction step.

<table><caption>Tab. 3. Hyperparameters exploration ranges for recurrent models by lead time and input forecast strategy. For the lead time (l) and sliding time windows (STW) given in numbers of time-steps (10 days)</caption><thead><tr><th>Recurrent model type</th><th>Lead time (l)</th><th>Number of neurons</th><th>STW ETRobs</th><th>STW ETRavg</th><th>STW Qobs</th><th>STW Qavg</th><th>STW Tobs</th><th>STW Tavg</th><th>STW RRobs</th><th>STW RRavg</th></tr></thead><tbody><tr><td rowspan="3">Rec</td><td>3</td><td>2 to 5</td><td>2 to 6</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 6</td><td>-</td></tr><tr><td>6</td><td>2 to 4</td><td>2 to 6</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 6</td><td>-</td></tr><tr><td>9</td><td>2 to 4</td><td>4 to 10</td><td>-</td><td>3 to 10</td><td>-</td><td>2 to 10</td><td>-</td><td>2 to 10</td><td>-</td></tr><tr><td rowspan="3">Rec_perfect</td><td>3</td><td>2 to 6</td><td>2 to 6</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 6</td><td>-</td></tr><tr><td>6</td><td>2 to 5</td><td>2 to 6</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 8</td><td>-</td><td>2 to 6</td><td>-</td></tr><tr><td>9</td><td>3 to 6</td><td>9 to 13</td><td>-</td><td>1 to 5</td><td>-</td><td>9 to 13</td><td>-</td><td>9 to 13</td><td>-</td></tr><tr><td rowspan="3">Rec_average</td><td>3</td><td>4 to 6</td><td>1 to 3</td><td>3</td><td>1 to 5</td><td>3</td><td>1 to 5</td><td>3</td><td>1 to 3</td><td>3</td></tr><tr><td>6</td><td>3 to 5</td><td>1 to 5</td><td>3</td><td>1 to 5</td><td>3</td><td>1 to 5</td><td>3</td><td>1 to 3</td><td>3</td></tr><tr><td>9</td><td>3 to 6</td><td>2 to 8</td><td>9</td><td>2 to 8</td><td>9</td><td>2 to 8</td><td>9</td><td>2 to 6</td><td>9</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 5. General architecture of the recurrent model.</figcaption></figure>

# 3. Results

## 3.1 Selected models

Based on the general architecture, the combined selection of input variables and model complexity led to the definition of final hyperparameters settings for each lead time (30, 60, and 90 days) and for each model type (Rec, Rec_perfect, and Rec_average). The optimal configurations for each forecast horizon, resulting from cross-validation selection based on the persistence criterion, are summarized in Table 4.

<table><caption>Tab. 4. Final hyperparameters configurations selected for each recurrent model type and forecast lead time. For the lead time (l) and sliding time windows (STW) given in numbers of time-steps (10 days)</caption><thead><tr><th>Recurrent model type</th><th>Lead time</th><th>Number of neurons</th><th>STW ETRobs</th><th>STR Q</th><th>STR T</th><th>STR RR</th><th>C<sub>p</sub></th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">Rec</td><td>3</td><td>5</td><td>3</td><td>4</td><td>6</td><td>2</td><td>0.41</td><td>0.31</td></tr><tr><td>6</td><td>2</td><td>2</td><td>3</td><td>6</td><td>4</td><td>0.64</td><td>0.23</td></tr><tr><td>9</td><td>4</td><td>5</td><td>9</td><td>2</td><td>8</td><td>0.72</td><td>0.37</td></tr><tr><td rowspan="3">Rec_perfect</td><td>3</td><td>5</td><td>4</td><td>2</td><td>3</td><td>6</td><td>0.72</td><td>0.73</td></tr><tr><td>6</td><td>4</td><td>5</td><td>6</td><td>4</td><td>6</td><td>0.87</td><td>0.67</td></tr><tr><td>9</td><td>3</td><td>10</td><td>5</td><td>9</td><td>12</td><td>0.91</td><td>0.58</td></tr><tr><td rowspan="3">Rec_average</td><td>3</td><td>4</td><td>6</td><td>6</td><td>8</td><td>4</td><td>0.38</td><td>0.42</td></tr><tr><td>6</td><td>3</td><td>5</td><td>5</td><td>6</td><td>4</td><td>0.65</td><td>0.42</td></tr><tr><td>9</td><td>6</td><td>13</td><td>17</td><td>17</td><td>13</td><td>0.69</td><td>0.26</td></tr></tbody></table>
The Rec_perfect model, benefiting from “ideal” future inputs, consistently shows the highest performance, highlighting the crucial role of future data in enhancing forecast accuracy. In contrast, the Rec model, using no future inputs, remains limited, with only slight improvement at longer horizons. The Rec_average model, relying on historical averages, provides only marginal gains at short lead times and significantly deteriorates beyond 30 days.

## 3.2 Performance of the models on the test set

After selecting the architectures for each forecast horizon, training is performed on the training set and stopped with the stop set. Generalization performance is evaluated on the test set. The test set results for the different models and forecast horizons are presented below. Table 5 shows the performance using $C_p$ and $R^2$ of the recurrent models Rec, Rec_perfect, and Rec_average. Fig. 7 provides a comparative evaluation of these models for 30, 60, and 90 days forecasting on the test set, with a graphical analysis of observed and predicted values.

<table><caption>Tab. 5. Performance of recurrent models at different lead times for the test set ($C_p$ and $R^2$ Metrics).</caption><thead><tr><th>Lead time</th><th>Recurrent model type</th><th>C<sub>p</sub></th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="4">30 days</td><td>Rec</td><td>0.054</td><td>0.31</td></tr><tr><td>Rec_perfect</td><td>0.8</td><td>0.84</td></tr><tr><td>Rec_average</td><td>0.09</td><td>0.45</td></tr><tr><td>Rec</td><td>-0.21</td><td>0.2</td></tr><tr><td rowspan="4">60 days</td><td>Rec_perfect</td><td>0.7</td><td>0.83</td></tr><tr><td>Rec_average</td><td>-0.45</td><td>0.07</td></tr><tr><td>Rec</td><td>0.1</td><td>0.39</td></tr><tr><td>Rec_perfect</td><td>0.66</td><td>0.78</td></tr><tr><td rowspan="2">90 days</td><td>Rec_average</td><td>-0.56</td><td>-0.04</td></tr><tr><td>Rec_average</td><td>-0.56</td><td>-0.04</td></tr></tbody></table>

The overall evaluation of the models across the full test set confirms the trends observed at the 30-day horizon. The Rec_perfect model consistently achieves high performance up to 90 days ($C_p \ge 0.66$, $R^2 \ge 0.78$), demonstrating its ability to produce reliable forecasts when provided with perfect future inputs. In contrast, both Rec and Rec_average show a marked decline in performance as the forecast horizon increases. The Rec_average model shows negative $C_p$ and $R^2$ at 90 days, indicating forecasts worse than the climatological mean.

Overall, the three model configurations can adequately reproduce the seasonal dynamics of groundwater levels, characterized by a gradual decline during the dry months (spring-summer) and a marked recovery during wetter periods (autumn-winter).

The 30 days lead time Rec_perfect configuration, which benefits from “perfect” forecasts of observed exogenous input variables up to the forecast horizon, delivers the closest match to the observed data, as expected. It shows an improved ability to capture rapid variations, particularly during groundwater level rises at the beginning and end of 2017.

The 30 days lead time Rec_average and 30 days lead time Rec configurations perform less well, with larger discrepancies from observations, especially during transition periods.

Furthermore, the results show that all three models with 30 and 60 days horizons are able to successfully anticipate the dry period of late 2016. Notably, the 30 days lead time Rec_perfect configuration is also capable of forecasting this dry episode up to a 90 days lead time. However, none of the model configurations can fully predict the 2017 drought, regardless of the forecast horizon (30, 60, or 90 days). This highlights a significant limitation in the models' ability to capture prolonged or more atypical low-flow events, possibly due to inadequate representation of summer recharge processes or groundwater system storage dynamics.

## 4. Discussion

This work explored the capacity of recurrent models to forecast groundwater levels under different assumptions regarding the availability of future meteorological inputs. The models were generally able to capture the main seasonal patterns, such as the gradual decline during dry periods and the recovery during wetter months. The 30 days lead time Rec_perfect configuration, which benefits from ideal knowledge of future inputs up to the forecast horizon, delivered the best match with observations, as expected, particularly during transitional phases.

However, given the study's main goal of predicting drought conditions -specifically low groundwater levels (e.g., in October 2016 and September 2017)- several limitations were identified. All models were able to anticipate the 2016 drought at 30- and 60-day horizons, and the Rec_perfect configuration extended this skill to 90 days. In contrast, none of the models successfully forecasted the 2017 drought, regardless of lead time.

This limitation could be attributed to an overreaction of the models to summer rainfall, which leads to premature and exaggerated increases in simulated groundwater levels. In practice, summer precipitation often fails to produce effective recharge due to high evapotranspiration rates and dry soil conditions. As a result, the models tend to overestimate system responsiveness, lacking the ability to capture delayed recharge and the inertia of groundwater dynamics.

Another limitation stems from the simplicity of the modeling framework. The architecture does not incorporate previous groundwater level measurements as additional inputs, which restricts its ability to produce forecasts from a known state. While this approach simplifies training and avoids feedback instability, it may reduce the model's effectiveness for medium and long-range forecasts.

Finally, despite the visual similarity of the predicted curves, performance metrics highlight substantial losses in accuracy when future inputs are missing or poorly represented. The Rec_average configuration, relying on historical climate means, consistently underperforms -especially beyond 30 days- demonstrating that long-term averages are not suitable substitutes for real or forecasted data in such a case. However, this is entirely consistent, since we are focusing in this study on a situation of extreme drought.

Models fed with “perfect” exogenous variables are the best regardless of the forecast horizon. This allows us to assess the importance of these variables, which not only provide the best forecasts in tests but also allow the model to catch the dynamics of the aquifer's behavior correctly during the training process. Improving the forecasting of these variables is therefore a major challenge for groundwater management.
<figure><img src="image_7.png" /><figcaption>Fig. 7. Observed vs. predicted piezometric levels using the three recurrent models. (a), (b), and (c) correspond respectively to 30-, 60-, and 90-days lead times.</figcaption></figure>

**5. Conclusion**

This study investigated the potential of simple recurrent neural network architectures to forecast groundwater levels at seasonal horizons (30, 60, and 90 days), with a particular focus on predicting drought episodes. The evaluated configurations differed in the availability and quality of future meteorological inputs, ranging from no future data to “perfect” forecasts.

Several key insights emerged from the results:

* Future meteorological data is critical: models provided with realistic or ideal future weather inputs clearly outperform those relying solely on historical data or long-term climatological means.

* Model simplicity limits forecast performance at extended lead times, especially when no mechanism is in place to incorporate previous groundwater level data as inputs.

* Drought events remain challenging to predict, particularly due to an over-sensitivity of the models to summer rainfall, which is often misinterpreted as effective recharge, despite high evapotranspiration and limited infiltration during dry periods.

* Historical averages are not an adequate proxy for future meteorological conditions in data-driven forecasting, especially beyond the short-term horizon and for extreme droughts.

While the predicted groundwater level trends align well with seasonal dynamics overall, the inability of the models to accurately reproduce low-level extremes—especially during critical drought periods—reveals structural limitations in both the data
representation and model design. In particular, the findings underscore the need for models to better capture the delayed and nonlinear nature of recharge processes, as well as the inertia of the groundwater system.

To address these challenges, several future directions can be explored:

* Enhancing rainfall input representation, by incorporating effective rainfall through preprocessing or dedicated model branches, could help filter out rainfall events with no hydrogeological impact and better reflect actual recharge contributions.

* Introducing observed piezometric level data as inputs in a feedforward architecture may improve temporal consistency and robustness, especially for short-term horizons like 30 days.

* Leveraging seasonal climate forecasts, even imperfect ones, offers a promising path toward operational improvements. The strong performance of the Rec_perfect model demonstrates the potential benefit of integrating weather forecast data into groundwater level prediction workflows.

These perspectives pave the way for the development of hybrid modeling approaches, combining meteorological forecasting, data processing techniques, and machine learning. Such models have the potential to deliver more reliable, generalizable, and actionable forecasts, particularly in support of proactive groundwater management under increasing pressure from climate variability and water demand.

## References

1. J. D. Mackay, C. R. Jackson, A. Brookshaw, A. A. Scaife, J. Cook and R.S. Ward, "Seasonal Forecasting of Groundwater Levels in Principal Aquifers of the United Kingdom", Journal of Hydrology, 530, 815-828 (2015).

2. I. N. Daliakopoulos, P. Coulibaly and I. K. Tsanis, "Groundwater level forecasting using artificial neural networks", Journal of Hydrology, 309 1-4, 229-240 (2005).

3. A. Wunsch, T. Liesch and S. Broda, "Forecasting Groundwater Levels Using Nonlinear Autoregressive Networks with Exogenous Input (NARX)", Journal of Hydrology, 567, 743-758 (2018).

4. A. Weisheimer and T. N. Palmer, "On the reliability of seasonal climate forecasts", Journal of The Royal Society Interface, 11(96):20131162 (2014).

5. A. Bukantis and V. Gytis, "Accuracy of monthly and seasonal forecasts generated for the territory of Lithuania using NOAA's Climate Forecast System version 2". Journal of Environmental Engineering and Landscape Management 29:337-45 (2021).

6. N. Wanders and E. F. Wood, "Improved Sub-Seasonal Meteorological Forecast Skill Using Weighted Multi-Model Ensemble Simulations", Environmental Research Letters 11(9):094007 (2016).

7. L. J. Slater, G. Villarini and A. A. Bradley, "Weighting of NMME Temperature and Precipitation Forecasts across Europe". Journal of Hydrology 552, 646-659 (2017).

8. M. Özger, A. K. Mishra and V. P. Singh, "Long Lead Time Drought Forecasting Using a Wavelet and Fuzzy Logic Combination Model: A Case Study in Texas" Journal of Hydrometeorology 13(1), 284-97 (2012)

9. P. Aghelpour, M. Babak, et S. M. Biazar. "Long-Term Monthly Average Temperature Forecasting in Some Climate Types of Iran, Using the Models SARIMA, SVR, and SVR-FA", Theoretical and Applied Climatology 138(3), 1471-1480 (2019).

10. A. Khastagir, I. Hossain, et A. H. M. Faisal Anwar, "Efficacy of Linear Multiple Regression and Artificial Neural Network for Long-Term Rainfall Forecasting in Western Australia", Meteorology and Atmospheric Physics 134(4):69 (2022).

11. Z. Wu, H. Yin, H. He and Y. Li, "Dynamic-LSTM Hybrid Models to Improve Seasonal Drought Predictions over China", Journal of Hydrology, 615A, 1287-1296 (2022).

12. S. Ju-Young, K. Kyu Rang and H. Jong-Chul, "Seasonal Forecasting of Daily Mean Air Temperatures Using a Coupled Global Climate Model and Machine Learning Algorithm for Field-Scale Agricultural Management", Agricultural and Forest Meteorology, 281, 1078-1089 (2020).

13. F. Habets, A. Boone, J. L. Champeaux, P. Etchevers, L. Franchistéguy, E. Leblois, E. Ledoux, P. Le Moigne, E. Martin, S. Morel, J. Noilhan, P. Quintana Segui, F. Rousset-Regimbeau and P. Viennot, "The SAFRAN-ISBA-MODCOU hydrometeorological model applied over France", Journal of Geophysical Research D: Atmospheres 113, D06113 18 (2008).

14. Bureau de Recherches Géologiques et Minières (BRGM) et Agence de l'Éau Seine-Normandie, "Fiche de caractérisation de la masse d'eau HG504", fiches de synthèse par masse d'eau, mars 2015.

15. K. Hornik, M. Stinchcombe, et H. White, "Multilayer feedforward networks are universal approximators", Neural Networks, vol. 2, no 5, p. 359-366, janv. 1989, doi: 10.1016/0893-6080(89)90020-8.

16. A. R. Barron, "Approximation and Estimation Bounds for Artificial Neural Networks", in COLT, M. K. Warmuth et L. G. Valiant, Éd., Oxford: Morgan Kaufmann, 1991, p. 243-249. doi: 10.1016/B978-1-55860-213-7.50025-0.

17. O. Nerrand, P. Roussel-Ragot, L. Personnaz, G. Dreyfus, et S. Marcos, "Neural Networks and Nonlinear Adaptive Filtering: Unifying Concepts and New Algorithms", Neural Computation, vol. 5, p. 165-199, mars 1993, doi: 10.1162/neco.1993.5.2.165.
References

18. S. Geman, E. Bienenstock, et R. Doursat, « Neural Networks and the Bias/Variance Dilemma », Neural Computation, vol. 4, no 1, p. 1-58, janv. 1992, doi: 10.1162/neco.1992.4.1.1.

19. J. Sjöberg et al., « Nonlinear black-box modeling in system identification: a unified overview », Automatica, vol. 31, no 12, p. 1691-1724, déc. 1995, doi: 10.1016/0005-1098(95)00120-8.

20. M. Stone, « Cross-Validatory Choice and Assessment of Statistical Predictions », Journal of the Royal Statistical Society: Series B (Methodological), vol. 36, no 2, p. 111-133, janv. 1974, doi: 10.1111/j.2517-6161.1974.tb00994.x.

21. T. Darras, A. Johannet, B. Vayssade, L. Kong-A-Siou, S. Pistre, « Ensemble model to enhance robustness of flash flood forecasting using an Artificial Neural Network: case-study on the Gardon Basin (south-eastern France) ». Boletin Geologico y Minero, 2018, 129 (3), pp.565-578. (10.21701/bolgeomin.129.3.007).

22. N. Akil, G. Artigue, S. Michaël, A. Johannet, et M. Vinches, « Uncertainty Estimation in Hydrogeological Forecasting with Neural Networks: Impact of Spatial Distribution of Rainfalls and Random Initialization of the Model », Water, vol. 13, p. 29, juin 2021, doi: 10.3390/w13121690.

23. Kitanidis, P.K.; Bras, R.L, "Real-time forecasting with a conceptual hydrologic model: 2. Applications and results", Water Resour. Res, 16, 1034–1044 (1980).

24. Nash, Je. and Sutcliffe, J. V, "River flow forecasting through conceptual models part I—A discussion of principles", J. Hydrol., 10(3), 282–290 (1970).