Paper deep dive
Pure and Physics-Guided Deep Learning Solutions for Spatio-Temporal Groundwater Level Prediction at Arbitrary Locations
Matteo Salis, Gabriele Sartor, Rosa Meo, Stefano Ferraris, Abdourrahmane M. Atto
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 97%
Last extracted: 3/31/2026, 1:33:22 AM
Summary
The paper introduces STAINet, an attention-based deep learning model designed for spatio-temporal groundwater level prediction at arbitrary locations. It incorporates physics-guided strategies—specifically inductive bias (STAINet-IB), learning bias (STAINet-ILB), and recharge zone information (STAINet-ILRB)—to improve model trustworthiness and generalization. The STAINet-ILB model demonstrated superior performance in predicting groundwater levels in the Piedmont region of Italy, effectively integrating sparse sensor data with dense weather information.
Entities (6)
Relation Signals (3)
STAINet → uses → Attention Mechanism
confidence 100% · we leveraged the attention mechanism
STAINet-ILB → trainedon → ERA5-land
confidence 95% · We adopted ERA5-land [32] as the source data for weather information.
STAINet-ILB → outperforms → STAINet
confidence 90% · The STAINet-ILB performed the best, achieving overwhelming test performances
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Groundwater represents a key element of the water cycle, yet it exhibits intricate and context-dependent relationships that make its modeling a challenging task. Theory-based models have been the cornerstone of scientific understanding. However, their computational demands, simplifying assumptions, and calibration requirements limit their use. In recent years, data-driven models have emerged as powerful alternatives. In particular, deep learning has proven to be a leading approach for its design flexibility and ability to learn complex relationships. We proposed an attention-based pure deep learning model, named STAINet, to predict weekly groundwater levels at an arbitrary and variable number of locations, leveraging both spatially sparse groundwater measurements and spatially dense weather information. Then, to enhance the model's trustworthiness and generalization ability, we considered different physics-guided strategies to inject the groundwater flow equation into the model. Firstly, in the STAINet-IB, by introducing an inductive bias, we also estimated the governing equation components. Then, by adopting a learning bias strategy, we proposed the STAINet-ILB, trained with additional loss terms adding supervision on the estimated equation components. Lastly, we developed the STAINet-ILRB, leveraging the groundwater body recharge zone information estimated by domain experts. The STAINet-ILB performed the best, achieving overwhelming test performances in a rollout setting (median MAPE 0.16%, KGE 0.58). Furthermore, it predicted sensible equation components, providing insights into the model's physical soundness. Physics-guided approaches represent a promising opportunity to enhance both the generalization ability and the trustworthiness, thereby paving the way to a new generation of disruptive hybrid deep learning Earth system models.
Tags
Links
- Source: https://arxiv.org/abs/2603.25779v1
- Canonical: https://arxiv.org/abs/2603.25779v1
Trouble viewing inline? Open PDF directly →
Full Text
89,340 characters extracted from source content.
Expand or collapse full text
Salis et al RESEARCH PAPER - PREPRINT Crossmark RECEIVED d Month y REVISED d Month y Pure and Physics-Guided Deep Learning Solutions for Spatio-Temporal Groundwater Level Prediction at Arbitrary Locations Matteo Salis 1,3,∗ , Gabriele Sartor 1 , Rosa Meo 1 , Stefano Ferraris 2 , Abdourrahmane M. Atto 3 1 Computer Science Department, University of Turin, Turin, Italy 2 Interuniversity Department of Regional and Urban Studies and Planning, Politecnico di Torino and University of Turin, Turin, Italy 3 LISTIC Laboratory, Universit ́e Savoie Mont Blanc, Annecy-le-vieux, France ∗ Corresponding author. E-mail: matteo.salis@unito.it or salis.matteo98@gmail.com (corresponding author) Keywords: Physics-Guided Machine Learning, Spatio-Temporal Data, Geospatial AI, Groundwater, Hydrology Abstract Groundwater represents a key element of the water cycle, yet it exhibits intricate and context-dependent relationships that make its modeling a challenging task. Theory-based models have been the cornerstone of scientific understanding. However, their computational demands, simplifying assumptions, and calibration requirements limit their use. In recent years, data-driven models have emerged as powerful alternatives. In particular, deep learning has proven to be a leading approach for its design flexibility and ability to learn complex relationships. We proposed an attention-based pure deep learning model, named STAINet, to predict weekly groundwater levels at an arbitrary and variable number of locations, leveraging both spatially sparse groundwater measurements and spatially dense weather information. Then, to enhance the model’s trustworthiness and generalization ability, we considered different physics-guided strategies to inject the groundwater flow equation into the model. Firstly, in the STAINet-IB, by introducing an inductive bias, we also estimated the governing equation components. Then, by adopting a learning bias strategy, we proposed the STAINet-ILB, trained with additional loss terms adding supervision on the estimated equation components. Lastly, we developed the STAINet-ILRB, leveraging the groundwater body recharge zone information estimated by domain experts. The STAINet-ILB performed the best, achieving overwhelming test performances in a roll- out setting (median MAPE 0.16%, KGE 0.58). Furthermore, it predicted sensible equation components, providing insights into the model’s physical soundness. Physics-guided ap- proaches represent a promising opportunity to enhance both the generalization ability and the trustworthiness, thereby paving the way to a new generation of disruptive hybrid deep learning Earth system models. 1 Introduction In recent years, deep learning models have achieved striking results in modeling many natural phe- nomena, including water resources [1, 2, 3, 4]. However, good predictions in terms of accuracy and precision are not enough. In some contexts, it is pivotal to have models that are not only able to replicate observed data but truly understand the generating process under study. Furthermore, black-box approaches can be challenging to employ when the goal is to advance scientific under- standing or support binding decision-making, given the lack of interpretation of what they have learned [5]. This poses major concerns in critical and high-stakes applications (e.g., natural hazards, healthcare, legal domain, etc.). One of the primary issues in studying hydrological and, more generally, environmental phenomena is the limited availability and quality of data. Indeed, environmental datasets are often incomplete, and labelled instances may fail to capture the full complexity of the underlying systems. Natural processes often exhibit intricate, non-stationary dynamics that evolve over time, influenced also by 1 arXiv:2603.25779v1 [cs.LG] 26 Mar 2026 Salis et al exogenous factors (e.g., climate change). Thus, inferring future trends from past observations might not be effective in some circumstances, because the process could change following unseen patterns – a situation well illustrated by the inductivist turkey 1 . For this reason, data-driven approaches, especially for environmental applications, might easily focus on spurious associations, making models unable to generalize well over the available spatio-temporal domain of the observations [5, 6]. Until now, theory-based models have been the milestone for our scientific knowledge and discov- ery. Briefly, theory-based approaches are based on cause-and-effect relationships, either empirically proven or theoretically demonstrated, starting from some postulates. Some examples are closed-form equations (e.g., continuity equations) or computational simulations of dynamical systems described by differential equations (e.g., Richards and Navier-Stokes equations). Although theory-based models remain invaluable for advancing our understanding of natural systems, they often rely on simplify- ing assumptions either because of incomplete knowledge of certain processes or to make the model computationally tractable. However, these assumptions are not always guaranteed in nature, and observations may diverge from model simulations [5, 7]. Furthermore, many internal parameters are usually required and thus estimated for the specific application, a procedure usually named model calibration. To calibrate a model, many observations are required; however, it is not al- ways feasible to obtain such data, either because of the expensive measurements or observational restrictions. This leads again to simplifying assumptions that may cause the model performance limitations [8, 5, 9, 10, 7]. Given that theory-based and data-driven approaches have their own advantages and disadvan- tages, a consistent number of researchers have proposed a mixed strategy to leverage the strengths of both, i.e., theory-guided data science [5, 6, 7], and specifically for physics and machine learning, the physics-guided (or informed) machine learning [11, 12, 13]. The idea is to inject prior scientific knowledge, a bias (or guidance), into a data-driven model, and to this aim, several approaches have been proposed, which differ in where to place the bias into the modeling pipeline. In general, three main strategies can be identified [12, 13]: observational bias, inductive bias, and learning bias. The observational bias strategy is conceptually the simplest one, as it simply assumes that the data considered truly and completely represent the underlying physical process. This could be achieved either by using accurate ground measurement, but as already pointed out, measurements could be scarce and costly, or by augmenting data with theory-based simulation [14]. The second strategy, the inductive bias, consists of designing the data-driven model in a way that implicitly embeds the domain knowledge. In the case of neural networks, it means developing specific layers that strictly enforce the desired relation. For example, as presented in [15] to enforce the advection-diffusion equation on sea-surface temperature data, and in [16] to find the viscosity solutions of two sets of Hamilton–Jacobi partial differential equations. Even though this strategy can constrain the model to strictly adhere to the domain prior, it requires ad hoc implementation, and usually, only simple and well-defined equations are injected [13]. The third strategy is the most general one, and it involves adding to the loss function between true and predicted data (L data ) some regularization terms (L eq ) related to the residuals of the prior gov- erning equation with respect to the model output. In general, the governing equation that describes the evolution of a phenomenon is expressed as a Partial Differential Equation (PDE). Within this framework, the model is trained to fit both the observed data and the physical constraints imposed by domain knowledge. This approach exploits the universal approximation capability of neural networks, enabling them to approximate PDE solutions leveraging automatic differentiation, and without requiring explicit domain discretization [11, 17]. Neural architectures designed for this pur- pose are commonly referred to in the literature as Physics-Informed Neural Networks (PINNs) [11]. The advantages of this last approach are: a) the neural network can be trained on the loss related to the prior equation (L eq ) with limited (or even any) labelled data; and b) if the differential equation contains some unknown parameters, these could be jointly estimated with the model’s own param- eters during the training process. However, given that the prior is injected as a regularization term into the loss, the model is not strictly, but softly bound to the domain knowledge, and thus could produce output not completely adhering to the prior equation. Although all three strategies have been demonstrated to be effective, the inductive and learning biases are the most adopted in recent years [12, 13, 6], the former because of the possibility of introducing hard constraints, the latter for its flexibility. In our application, we focused on groundwater, one of the major freshwater components in the water-cycle, whose accurate accounting is critical for water policy planning, especially in the presence 1 The inductive turkey (also known as the turkey illusion) exemplifies the problem of induction, originally described by Bertrand Russell in The Problems of Philosophy. It tells of a turkey that, being fed daily, infers that this pattern will persist indefinitely—until, on the day before Thanksgiving, it is slaughtered, exposing the weakness of its assumption. 2 Salis et al of climate change. In more detail, this work addressed two main objectives. Firstly, we aimed to model the weekly groundwater level in Piedmont (Italy), evaluating different strategies for the injection of the ground- water flow equation. Additionally, given our interest in developing a general model for the whole Area of Interest (ROI), we sought to build models able to predict at any desired (i.e., arbitrary) location all over the ROI, leveraging an autoregressive and an exogenous component, which consist of the spatially sparse groundwater measurements derived from in situ sensors (piezometers) and the spatially dense, or distributed, weather information (weather video). To the best of our knowledge, this is the first study to test different physics-guided deep learning approaches in a real-world sce- nario, rather than on simulated groundwater data. 1.1 Groundwater principles Our study addressed the shallow groundwater bodies, named also phreatic aquifers or unconfined aquifers, which are portions of the terrain saturated with water and directly influenced by the superficial activity, and whose upper surface, denoted as the water table, is at atmospheric pressure. To quantify groundwater resources, the groundwater level (GWL) is frequently adopted, and it represents the distance between the reference datum and the water table. Groundwater dynamics are usually described by using Darcy’s law, which states that the specific discharge q x in the x-direction is proportional to the gradient of the total hydraulic head h multiplied by the saturated hydraulic conductivity K x . More specifically, the specific discharge q x is the volume rate of flow per unit of area, and the hydraulic head is the potential mechanical energy per unit weight of fluid. The hydraulic head is the sum of the gravitational head, given by the elevation of the point above an arbitrary horizontal datum, and the pressure head, which is the pressure of the fluid at the point divided by the specific weight of the fluid. In the case of an unconfined aquifer, whose upper surface is at atmospheric pressure, the hydraulic head is equivalent to the water table height (i.e., groundwater level). Therefore, it is possible to write Darcy’s law as in Equation 1: q x =−K x ∂h ∂x (1) Combining Darcy’s law with the mass conservation principle, it is possible to derive a general PDE, a form of diffusion equation, for the saturated groundwater flow in the three spatial dimensions (x,y,z), as in Equation (2). The term S s represents the specific storage 2 [L -1 ], that is the volume of water entering/leaving from a unit volume per unit variation in hydraulic head, and K x , K y , K z are the conductivity in the specific directions, all expressed in [LT -1 ]. ∂h ∂t S s = K x ∂ 2 h ∂x 2 + K y ∂ 2 h ∂y 2 + K z ∂ 2 h ∂y 2 (2) However, attention is typically directed toward horizontal flow along the aquifer plane of coordi- nates (x,y). Considering a simplified formulation and defining the transmissivity T = H·K [L 2 T -1 ] and the storage coefficient S = H · S s [1], in which H [L] is the saturated thickness of the aquifer and K is the horizontal conductivity, it is possible to define the groundwater flow equation for the horizontal flow (hereafter groundwater flow equation) as in Equation (3): ∂h ∂t S = T ∂ 2 h ∂x 2 + ∂ 2 h ∂y 2 (3) We reformulated Equation (3), moving S to the right-hand side obtaining Equation (4). The resulting quantity T S is related to the ease with which water moves through the soil, and it is expressed in [L 2 T -1 ]. ∂h ∂t = T S ∂ 2 h ∂x 2 + ∂ 2 h ∂y 2 (4) 2 Related Works An increasing number of studies have explored the implementation of physics-guided approaches across a wide range of domains, including healthcare [18, 19], computer vision [20], geoscience [15, 21, 22, 23], weather [24, 25], air-pollution [26], and hydrology [27, 22, 28, 29]. 2 As common in hydrology, to describe a unit of measurement, we adopted the terminology of the fundamental physical dimensions, thus length [L] and time [T] – in our case meters [m] and weeks [w] 3 Salis et al More in detail, the work presented in [15] aimed to develop a deep learning model to predict sea surface temperature in an autoregressive fashion, leveraging the information of the advection- diffusion equation and adopting both inductive and learning bias strategies. The authors developed a specific architecture made of two branches. The first, made of a U-Net-like structure, explicitly estimates the motion field; the second enforces the motion field on the last available temperature image through a mathematically consistent warping scheme. The warping scheme was developed to obtain the demonstrated discretized solution using a radial basis function kernel (inductive bias). The relevant point in this work is that there is no direct supervision on the motion vector field, but only a weak supervision through some regularization terms added to the final loss (learning bias). With this model, the authors achieved better performance than a numerical assimilation model and other deep learning implementations, proving the effectiveness of their solution. Concerning hydrology, the inductive bias strategy was adopted in [27] to develop a physics- guided neural network with the aim of modeling runoff across the conterminous United States. In particular, they proposed an architecture comprising two consecutive blocks. The first is composed of a restructured version of an RNN layer, named P-RNN, which embeds the spatially lumped hydrological model described by two discrete state-space representations. Then the second block, made of usual 1D convolutional layers, takes the original input along with the P-RNN output to model the residual component with respect to the physics. The proposed physics-guided model achieved better performance compared to other deep learning implementations. A different strategy was followed by authors in [28, 29]. Indeed, they adopted the PINNs approach (learning bias) to estimate groundwater flow simulated data. In the first study [28], a multilayer- perceptron (MLP) was trained to estimate the solution of the groundwater flow equation in one and two-dimensional cases in the presence of pumping wells, by adding to the loss a residual term of the PDE and boundary condition terms. The trained PINN was revealed to be more accurate than the finite difference method adopted as a comparison. In the second study [29], the authors aimed to estimate the phreatic surface (i.e., water-table) and the piezometric heads in a vertical cross-section of a homogeneous and isotropic aquifer and a heterogeneous and anisotropic one. To this aim, they proposed a physics-informed neural network made of two sub-networks, each implemented as an MLP. The first is responsible for estimating the piezometric head value along the vertical cross- section over time, while the second network is designed to predict the phreatic surface. In that study, MODFLOW [30] simulations were performed and considered as the true data in the training process, jointly with the residuals of the groundwater flow equation. Even though the architecture appears to be intricate, the authors achieved better performance than a pure deep learning architecture without the physics loss. Furthermore, they highlighted the ability of PINNs to work in mesh-free domains (i.e., without discretization) and that incorporating physical constraints can dramatically reduce the number of required observations for training. This shows the PINNs’ potential to generalize well also in data-scarce environments. Even if these mentioned studies reveal the potential of a physics-guided deep learning approach, our application context is different. Indeed, in [15], the adopted target (sea surface temperature) is dense in space and time. This contrasts with our case, where groundwater level observations are spatially sparse, and predictions may be required at locations different from those of the input data. On the other hand, even if the PINN approach is flexible, it is usually implemented for theoretical analyses, and thus models are frequently developed to take as input only the coordinates of the system (spatial, temporal, or both), and less frequently exogenous, and real, inputs are considered. Furthermore, feeding multi-modal data (e.g., video and tabular data) to the models in this framework is non-immediate, and a careful design is needed. In our study, to inject the groundwater flow equation, we implemented and tested two of the above mentioned strategies. Firstly, we introduced an inductive bias into the pure deep learning model by making the model estimate the different components of the governing equation. Afterwards, we added learning biases into the loss function to further condition the model. Concerning our other objective of producing models that can predict at arbitrary spatial and temporal locations independently of the locations at which the input features are observed, we leveraged the attention mechanism [31]. Indeed, attention maps values indexed by keys to new values indexed by queries (in a database fashion), which offers a straightforward solution to leverage input data indexed by spatio-temporal coordinates (keys) to generate predictions at an arbitrary number of new locations (queries). Other works, as [3, 25], leverage spatially sparse data as well; however, in the first case [3] authors interpolated data into a grid, while in the second case [25], the authors adopted an MLP to perform the upsampling, and thus the number of input and output locations must be fixed in advance. It is worth noting that developing models that automatically predict at different locations than 4 Salis et al the input ones is in contrast with traditional approaches, which require input features to be co- located with the target variable and then use post hoc interpolation methods to generate predictions in new locations. 3 Data In this study, we defined our ROI as a squared region with longitudes ranging from 6.6267° to 8.0292° and latitudes from 44.3108° to 45.2650°, of approximately 16700 km 2 (see Figure 1). In this area, we retrieved 28 piezometers’ weekly time series, measuring the groundwater level from 2001-01-01 to 2023-12-31. These sensors are part of the Regional Environmental Agency of Piedmont (ARPA Piemonte in Italian) network 3 . The time series of sensors in Cuneo, Racconigi, Scalenghe, and Vottignasco are reported in Figure 2; all the other time series, along with their means, standard deviations, and percentages of missing values, are shown in Section A in Figure 20 and Table 4. The right side of Figure 1 represents the location of all the sensors with a colour scale indicating their mean over time. A visible aspect of these time series is that they show very different dynamics and missing periods. This is a complex issue to handle because it suggests the presence of different processes involved in the groundwater level data. In particular, some exogenous phenomena, like abstraction for irrigation or the response to precipitation, might act differently in different locations. Given the high percentage of missing values among all sensors (see Table 4 in Section A) with an average value of around 23%, we decided not to discard any sensors but simply fill the missing data with the last available data for every sensor. These filled data were then masked out and thus not considered in the loss computation during training. We adopted ERA5-land [32] as the source data for weather information. ERA5-land is a reanalysis dataset providing different environmental variables at a spatial resolution of 0.1° and an hourly temporal resolution. In particular, we downloaded total precipitation, 2 m height temperature, potential evaporation, and snowmelt, and performed a weekly average aggregation. We considered data up to 2021-12-31 for training (around 90% of all instances), leaving data in 2022 and 2023 for testing (roughly 10%). Furthermore, we normalized all data by computing z-scores z = F−μ F σ F for all features F , where all means μ F and standard deviations σ F were computed on the training set. Figure 1: Piedmont region with the ROI box and the selected piezometers. 4 Methods Our target variable is the hydraulic head h(x,y,t) indexed by time t, and spatial coordinates x,y. Given that we considered the unconfined aquifer, h corresponds to the groundwater level (GWL) 3 Free access is provided at https://shorturl.at/F3KzQ. 5 Salis et al (a) (b) (c) (d) Figure 2: Time series of the sensors in a) Cuneo, b) Racconigi, c) Scalenghe, and d) Vottignasco. measured by piezometers. The ROI, or spatial domain, is Ω ⊂R 2 [ROI] delimited by the bounding box defined by (x min ,y min ), (x max ,y max ). The piezometers measure h at M = 28 points inside the ROI, we can thus define the spatial measuring domain M as the set of sensors’ locations M = (x m ,y m )|m∈ [1;M ] and M⊂ Ω. We aimed to develop a model that can predict an arbitrary number P of points P inside the ROI at time t ∗ , more precisely P = (x p ,y p )|p ∈ [1;P ]. Our target can thus be defined as h P,t ∗ = h(x,y,t)|(x,y)∈P,t = t ∗ , note that P ⊂ Ω. As input data, we leveraged an autoregressive component H P,t ∗ −1 made of T lagged values of the target, i.e., H P,t ∗ −1 =h P,t |t∈ [t ∗ −T ;t ∗ − 1]. Furthermore, we used spatially dense weather data all over the ROI in the form of a video, from the time of the former T lag up to the time t ∗ of the prediction. We defined the video as V t ∗ =v(x,y,c,t)∈R|(x,y)∈H×W,c∈ [1;C],t∈ [t ∗ −T ;t ∗ ], where H,W are the horizontal and vertical centroids of each pixel of the video that covers the entire spatial domain Ω, and C are the channels, i.e. total precipitation, 2 m height temperature, potential evaporation, and snowmelt. Formally, h P,t ∗ = f (H P,t ∗ −1 , V t ∗ ) + ε P,t ∗ is the considered theoretical equation generating GWL data, where ε P,t ∗ is the irreducible error. The relation f can be approximated by a model f θ parametrized by θ, in our case, a neural network. 6 Salis et al 4.1 Deep Learning Approach 4.1.1 Multi-Head Attention Mechanism To develop a deep learning model that leverages input data indexed by spatio-temporal coordinates and can generate predictions at an arbitrary and vari- able number of locations, we adopted the Multi-Head Attention (MHA) mechanism [31]. The MHA is based on the scaled dot-product attention that computes attention weights by taking the dot product between query (Q) and key (K) vectors, scaled by the inverse square root of the key dimension (d k ) to maintain numerical stability. These weights are normalized using a softmax function and applied to the corresponding value vectors (V ) to produce a weighted sum representing the attention output (Equation 5). When the queries, keys, and values are all derived from the same input sequence, the mechanism is referred to as self-attention. This mechanism enables efficient computation of dependencies within a sequence by emphasizing relevant contextual information while mitigating gradient instability in high-dimensional spaces. Attention(Q,K,V ) = softmax QK T √ d k V(5) The Multi-Head formulation consists of performing n parallel and independent scaled dot-product attention (heads) on different linear transformations of the input query Q i = QW Q i , value V i = V W V i , and key K i = KW K i for i ∈ [1;n]. All heads’ outputs are then concatenated, and a linear transformation is performed to produce the final output, as in Equation (6): Multi-Head Attention(Q,K,V ) = Concat(head 1 ,..., head n )W O (6) Where W Q i , W V i , W K i and W O are learnable weight matrices, and the generic head i is equal to Attention(Q i ,K i ,V i ). This structure enables the model’s capacity to extract different, but complementary, meaningful information from the input sequence and thus learning more complex relationships. The MHA layer is usually adopted in blocks, in which it is followed by feed-forward layers as in the Transformer architecture [31, 33]. For our purpose, and inspired by [33], we defined an MHA block as depicted in Figure 3a that consists of an MHA layer followed by an MLP made of two fully connected layers. We used the In- stance Normalization [34] as the normalization procedure (Norm layer in Figure 3a), and LeakyReLU as the activation function after each layer. Two residual connections are then employed to ease the gradient flow. 4.1.2 Embedding Layers The values, keys, and queries were obtained through dedicated embedding layers. Conceptually, in our context, the keys or queries refer to the spatio-temporal coordinates of the input data and the prediction points, respectively. To create a more informative spatio-temporal coordinates embeddings (ST-CEmb), we augmented the coordinates (longitude and latitude) with the elevation, along with the sine and cosine signals of the day of the year (doy) as the temporal information – the functions sin( 2π doy 366 ) and cos( 2π doy 366 ) respectively. The usage of sine and cosine signals for the embedding creation has been adopted in the original Transformer architecture [31] to inform about the position of each token in the sequence, but also in environmental applications to ease the seasonality extraction [35, 36]. The augmented spatio-temporal coordinates of each point, either in input or to be predicted, were fed to a linear layer that gives in output the spatio-temporal coordinates embeddings – the keys, or the queries, respectively. Concerning the computation of the spatio-temporal value embeddings (ST-VEmb), for the au- toregressive component, we adopted a linear layer that takes as input the groundwater level value concatenated with the augmented spatio-temporal coordinates described before. Differently, for the value embedding of the weather data, we exploited the spatially dense nature of this input by feeding the weather features, concatenated with the augmented spatio-temporal coordinates, into a convo- lutional layer with a 3×3 kernel, followed by an average pooling layer. Consequently, the resulting ST-VEmb inherently captures the spatial structure of the meteorological variables. 4.1.3 Deep Learning Architecture Upon the MHA block, we defined a) the Spatio-Temporal Self- Attention (STSA, in Figure 3b) module made of one MHA block adopting self-attention; and b) the Spatio-Temporal Attention-based Interpolation (STAI, in Figure 3b) module composed of one MHA block in which the keys are the spatio-temporal coordinates embeddings of the input points, the values are the hidden representation of the value embeddings, and the queries are the spatio- temporal coordinates embeddings of the prediction points – given that the keys and queries could be of different lengths, we dropped the first residual connection in the MHA block. 7 Salis et al (a) (b) (c) Figure 3: a) Multi-Head Attention block adopted in the models. b) STSA Module c) STAI Module that takes as values the input points hidden representation of the value embeddings, as keys the input points spatio-temporal coordinates embeddings (ST-CEmb), and as queries the prediction points spatio-temporal coordinates embeddings (ST-CEmb). Figure 4: STAINet architecture. To further condition the model on the spatio-temporal location of the prediction points, we adopted the Feature-wise Linear Modulation (FiLM) conditioning [37]. The FiLM layer performs a conditional feature normalization by scaling and shifting every feature channel. In particular, the scale and shift factors are computed by a linear layer from the conditioning information – in our case, the augmented spatio-temporal coordinates of the prediction points. This mechanism enables adaptation of neural representations to the conditioning information, facilitating the learning of complex context-dependent relations and improving the generalization ability of the model. We named our pure deep learning model Spatio-Temporal Attention-based Interpolation neural 8 Salis et al Network (STAINet), and its architecture is depicted in Figure 4. We decided to process the au- toregressive component and the weather one, with two parallel branches. In particular, we designed an Autoregressive Branch that begins with one STSA module. During training, the MHA layer adopted in this first module employs dropout to randomly omit part of the input key–value pairs, thereby enhancing the model’s generalization capability and preventing it from focusing on specific input locations [3]. Subsequently, to extract the relevant information from the input points, an STAI module projects the hidden representation of the spatio-temporal value embedding (i.e., the output of the previous STSA module) onto the spatio-temporal prediction query, followed by an additional STSA module. The FiLM conditioning is adopted both after the STAI module and after the last STSA module to force the hidden representations to capture context-dependent relationships and facilitate the model in differentiating among prediction points. This configuration was found to be the optimal choice for the Autoregressive Branch, balancing models’ complexity with the ability of learning complex dependencies from in situ spatially sparse measurements. To handle the exogenous weather information, we defined the Weather Branch, which is made only of one STAI module, whose outputs are conditioned by the FiLM layer. The inclusion of STSA modules in this branch did not prove effective; on the contrary, it made the training process more difficult and reduced the models’ generalization capability. This behaviour is likely attributable to the use of the time-distributed convolutional embedding layer for constructing the ST-VEmb, which already encapsulates the relevant spatial information. For the pure deep learning model, the two branches were concatenated, and the resulting hidden representation was passed through four stacked STSA modules. The output is then processed by a last linear layer that produces the predictions for all the desired points. Due to differences in the time series across sensors, we adopted the Mean Absolute Percentage Error (MAPE) as the data loss term (L data ) since it proved more effective than other non-relative loss functions. Furthermore, we introduced an orthogonality regularization term L ortho to enforce orthogonality on the weights of the embedding layers. This constraint encouraged the embeddings to capture non-redundant features [38, 39]. 4.2 Physics-Guided Approach Considering the PDE in Equation (4), applying the Euler integration method, we then expressed h at a time t ∗ (h t ∗ ) as a function of its ∆t lag (h t ∗ −∆t ), as in Equation (7). h t ∗ = h t ∗ −∆t + ∆t T S ∂ 2 h ∂x 2 + ∂ 2 h ∂y 2 (7) To simplify the notation, we defined the diffusion parameter D = T S expressed in [L 2 T -1 ], and ∆ GW t ∗ = ∆t h D ∂ 2 h ∂x 2 + ∂ 2 h ∂y 2 i expressed in [L]. We further add a residual (sink/source) term R t ∗ expressed in [L] related to all exogenous factors occurring in the interval ∆t, which can be either anthropogenic (e.g., water abstraction for irrigation) or natural (e.g., rainfall or snowmelt recharge). We thus obtained Equation 8: h t ∗ = h t ∗ −1 + ∆ GW t ∗ +R t ∗ (8) Hereafter, we refer to h t ∗ −1 as the autoregressive component and to ∆ GW t ∗ as the diffusion displacement component. In our setting, we considered ∆t equal to one week. 4.2.1 Inductive Bias Strategy As an initial approach to incorporate physical knowledge into the deep learning model, we modified STAINet to explicitly output the three components defined earlier, computing the final output as the sum of these components. To this end, we developed the PSTAINet- IB (Physics-guided STAINet with Inductive Bias), which is depicted in Figure 6. The initial part of the architecture follows the same structure as the STAINet; however, PSTAINet- IB does not concatenate the two branches. In particular, the R t ∗ and h P,t ∗ −1 terms are estimated by two linear layers. Differently, since the diffusion component ∆ GW t ∗ is more complex, we employed a more advanced neural structure for its estimation. Firstly, theD Block estimates the parameterD, taking as input the spatial coordinates and eleva- tion of the prediction points. It is implemented as an MLP comprising three linear layers (Figure 5a), followed by Instance Normalization and LeakyReLU activation. We designed the D Block under the assumption thatD is a time-invariant property of the soil, and thus can be estimated based on spatial coordinates. Although simple in structure, this block provides hydrologically consistent estimations of D while keeping the number of trainable parameters low. We then developed the ∆ GW module 9 Salis et al (a) (b) Figure 5: a) D Module in which Act & Norm stand for activation and normalization respectively; b) ∆ GW Module. Figure 6: PSTAINet-IB architecture. (Figure 5b), consisting of a stack of two STSA modules, to estimate the divergence term ∂ 2 h ∂x 2 + ∂ 2 h ∂y 2 . In the ∆ GW module, we adopted STSA modules in place of MLP to enable the estimation of more complex relationships among prediction points. The diffusion component is then computed as the multiplication of the estimated ˆ D and divergence term (∆t is omitted as we set it equal to one). During the training of the model PSTAINet-IB, the three terms were not directly supervised; only their sum, corresponding to the final output of the model, was used in the calculation of the data loss function L data . 4.2.2 Learning Bias Strategy We then tested a second strategy consisting of adding learning biases to the PSTAINet-IB. Firstly, we supervised the estimated autoregressive component using the MAPE between ˆ h P,t ∗ −1 and the true data h P,t ∗ −1 , to make the lag estimation coherent with the available data, we named this loss term L coh . Subsequently, we computed the Mean Squared Error (MSE) of the residuals between the estimated diffusion displacement and its discretized version computed on the estimated lag – we named the 10 Salis et al loss built on these residuals L diff . In detail, we applied the finite difference method with a central difference approximation [40] to compute the second-order spatial derivatives of ˆ h P,t ∗ −1 . To compute this discretized approximation, a spatially dense estimated lag is needed, i.e., a gridded map. It is worth noting that this loss can be computed on arbitrary prediction points, named control points in the PINN literature, and does not require true data for supervision. We further introduced some regularization terms on the ℓ 1 and ℓ 2 norms of the ˆ R to constrain its mean value toward zero and mitigate high-magnitude responses, termed L ||R|| 1 and L ||R|| 2 , respec- tively. We named the model with these learning biases PSTAINet-ILB (Physics-guided STAINet with Inductive and Learning Biases). Finally, we leverage the domain knowledge of recharge zones, which are the zones where the groundwater body recharges due to precipitation and snowmelt. In our ROI, ARPA experts defined the recharge zones for the groundwater body in an area at the base of the mountain 4 . We considered a buffer of 0.05° and defined a loss termL RCH that penalizes positive values ofR outside the recharge zones. Formally, definingR + as the positive value ofR, RCH as the set of points inside the recharge zones, and its complement RCH c , L RCH = 1 P 1 RCH c (R + ) P 1 RCH c (R + )R. We named the model trained with this further regularization PSTAINet-ILRB (Physics-guided STAINet with Inductive, Learning, and Recharge zone Biases). The final loss is then the summation of all terms previously described (if adopted by the considered model), as in Equation (9): L =L data + α ortho L ortho + α coh L coh + α diff L diff + α ||R|| 1 L ||R|| 1 + α ||R|| 2 L ||R|| 2 + α RCH L RCH (9) Table 1: Hyperparameters. Modelα ortho α coh α diff α ||R|| 1 α ||R|| 2 α RCH STAINet10 −3 ----- PSTAINet-IB10 −3 ----- PSTAINet-ILB10 −3 12.5· 10 −2 5· 10 −4 10 −4 - PSTAINet-ILRB10 −3 12.5· 10 −2 5· 10 −4 10 −4 5· 10 −4 5 Experiments We supplied the models with information from the four weeks preceding each prediction (lagsT = 4), as this configuration yielded the best performance while avoiding an excessive amount of uninforma- tive data. For all four models STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB, we set the embedding dimension equal to 32 and the number of heads for each MHA layer to 8. Under these settings, STAINet contained a total of 62081 trainable parameters, while PSTAINet-IB, PSTAINet- ILB, and PSTAINet-ILRB each comprise 49666 trainable parameters. All models were trained using AdamW optimizer [41] with a learning rate of 2.5· 10 −4 decreased to 10 −4 at epoch 400, and a weight decay equal to 2.5· 10 −3 . The total number of epochs to achieve a relevant optimum was found to be 1000 for STAINet, 850 for PSTAINet-IB and PSTAINet-ILRB, while 600 for PSTAINet-ILB. The dropout rate in the initial Spatio-Temporal Self-Attention block of the Autoregressive Module was set to 0.15 for all configurations. The weights assigned to the individual loss components are reported in Table 1. All the mentioned hyperparameters were found by manual grid search. We compared the model performances considering the median over all sensors of NBIAS, RMSE, MAPE, Nash–Sutcliffe Efficiency (NSE), and Kling–Gupta Efficiency (KGE), described in Equations 10, 11, 12, 13, and 14, computed between true data h and the prediction ˆ h for every sensor over the test set of length N = 104 (weeks), and with σ and μ representing the variance and mean respectively. To evaluate further the model’s generalization ability, we tested the models feeding both the true lags as input and in a rollout setting using the predictions at the previous time step as the lagged input for the whole test, i.e. ˆ h P,t ∗ = ˆ f θ (H P,t ∗ −1 , V t ∗ ) and ˆ h P,t ∗ = ˆ f θ ( ˆ H P,t ∗ −1 , V t ∗ ), respectively. Thus, in the last case, models see only true groundwater level data at the beginning of the test set (2022-01-01), and then they iterate their predictions to produce the whole test set’s predictions. Since the model can generate predictions at an arbitrary number of points, it enables the produc- tion of groundwater level maps at any desired spatial resolution (as shown in the output of Figure 4). 4 Available here: https://shorturl.at/qTheu. 11 Salis et al Therefore, we produced dense prediction maps at a spatial resolution of 1.5km per pixel. Specifically, prediction points were placed approximately every 1.5km, and all points were jointly predicted by the models. This allowed us to assess both the coherence of the predictions and the behaviour of the predicted equation components throughout the ROI area 5 . Finally, the best-performing model was used to predict the full time series (starting from 2001- 01-01) for all sensors, assessing its capability to replicate the observed data. In this last case, we produced rollout predictions with a forecast horizon of 26 time steps, i.e., feeding true data every 26 time steps, and iterating the predictions in between. RMSE = v u u t 1 N N X i ( ˆ h i − h i ) 2 (10) NBIAS = 1 N P N i ( ˆ h i − h i ) h max − h min (11) MAPE = 100 1 N N X i | ˆ h i − h i | h i (12) NSE = 1− P N i ( ˆ h i − h i ) 2 P N i (h i −h) 2 (13) KGE = 1− p (ρ− 1) 2 + (α− 1) 2 + (β− 1) 2 α = σ ˆ h σ h ;β = μ ˆ h μ h (14) 6 Results Table 2 reports the median evaluation metrics over all sensors computed on models’ predictions feeding real data and in the rollout setting. In both cases, the best-performing model was PSTAINet- ILB. In general, all models with the physics prior performed better than the pure deep learning approach (STAINet), suggesting that both the inductive bias and the learning bias strategies provided useful information to the models. On most of the metrics, PSTAINet-ILRB performed better than the PSTAINet-IB, but worst than the PSTAINet-ILB. This might suggest that the recharge zones prior was too restrictive, limiting the model’s representation ability. The metrics computed on every sensor and every model in the rollout setting are reported in Section B, specifically in Table 5, Table 6, Table 7, Table 8. Figure 7 depicts the predictions on the test set, feeding true data and iterating the predictions (rollout setting), respectively, for the sensors in Cuneo, Racconigi, Scalenghe, and Vottignasco. In Section B, we reported the predictions over the test set for all the sensors, in the rollout setting (Figure 21). From Figure 7, the good generalization ability of the PSTAINet-ILB model is visible, especially in contexts of unprecedented drop, like in the summers of 2022 and 2023. PSTAINet-ILRB also predicted the decreases in the water level among all sensors accurately, but in some cases (like in Scalenghe, Figure 7f) the drops were not accurately predicted, and in other cases (like Racconigi, Figure 7d) a premature recharge was predicted. Even if differences between feeding true data and the rollout setting are visible (Tables 2 and from Figure 7, the PSTAINet-ILB model maintained almost the same performance. This further demon- strates the good generalization ability of the PSTAINet-ILB model. Figure 8 shows the map predictions at a resolution of 1.5km/pixel for the week starting on 2022- 07-17. All models predicted a map in line with the expert domain expectations (i.e., with higher groundwater levels in the mountain region). However, the PSTAINet-ILB and PSTAINet-ILRB predictions appear to be less definite on the north-west mountains, but they better identified the rivers and valleys in the southeastern part. Notwithstanding, a formal evaluation of these results is difficult because of the lack of true measured data all over the ROI, especially in the mountains – a universal problem in hydrology. 5 The code is available here https://github.com/Matteo-Salis/physics-guided-gwl along with GIFs representing the prediction maps evolution through time 12 Salis et al (a) (b) (c) (d) (e) (f) (g) (h) Figure 7: STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB predictions on the test set feeding true data in a) , c) , e) , g) and in the rollout setting b) , d) , f ) , h) in Cuneo, Racconigi, Scalenghe, and Vottignasco. 13 Salis et al Table 2: Median performance metrics on Test set feeding true data (True Data) and in the rollout setting (Rollout). STAINetPSTAINet-IBPSTAINet-ILBPSTAINet-ILRB True Data Rollout True Data Rollout True Data Rollout True Data Rollout NBIAS0.10420.17130.05940.0862-0.00930.01240.01770.0685 RMSE [m]0.83511.02920.62960.73360.50580.58510.65210.6783 MAPE [%]0.28340.28680.23190.26060.14790.16040.19000.1953 NSE0.56120.28840.70590.62690.80250.78760.79860.7024 KGE0.25830.22360.42530.51960.59170.58000.50180.4074 Figure 8: STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB groundwater level predic- tion maps at 1.5km/pixel for 2022-07-17. 7 Discussion An extensive comparison with existing studies is challenging due to the specificity of our application and the scarcity of research employing physics-guided strategies on real groundwater level measure- ment data. However, a reference work exists [42], in which authors trained sensor-specific pure data-driven neural networks on a subset of groundwater sensors considered in our study. Although our proposed models are not sensor-specific, we achieved comparable performance, considering the PSTAINet-ILB median NBIAS −0.01 versus −0.05 of [42], NSE 0.80 versus 0.90, and KGE 0.59 versus 0.43. Our performance are consistent also with those reported in other studies that adopt data driven models to estimate groundwater levels in different geographical regions, such as [43] in Yangpyeong riverside area (South Korea) and [44] in the Upper Rhine Graben area (central Europe), which achieved NSE values of nearly 0.8 and 0.5, respectively. In the following, a more detailed analysis of the results of our proposed models is presented. 7.1 Loss Analyses Figures 9a and 9b, which report the values of L data on the training and test sets, indicate that PSTAINet-ILB and PSTAINet-ILRB exhibit superior generalization ability. Specifically, although both models obtained higher training L data values, they achieved lower test errors, demonstrating better out-of-sample performance. Furthermore, even if PSTAINet-ILB has been trained for fewer epochs (600) it achieved the best performance, further proving the utility of the adoption of both inductive and learning bias strategies. Nevertheless, the introduction of learning biases substantially increased the time required for completing a training epoch, especially due to gradient computations, which was on average 4.616s for STAINet, 4.157s for PSTAINet-IB, 18.286s for PSTAINet-ILB, and 18.697s for PSTAINet-ILRB. Figures 9c, 9d, and 9e report the values ofL coh ,L diff , andL RCH on the training set, respectively. Both L coh and L RCH exhibit a monotonic decrease over the course of training, whereas L diff shows an initial decline followed by a subsequent rebound for both models. This behaviour may stem from interactions or conflicts among the different loss terms, but it may also be attributed to noise, given the very small magnitudes reached by L diff (see the y-axis in Figure 9d). A more in-depth analysis of these loss dynamics is left for future work. 14 Salis et al (a)(b) (c)(d) (e) Figure 9: Loss plots a) training L data , b) test L data , c) training L coh , d) training L diff , e) training L RCH . 7.2 Error Analyses 7.2.1 Missing Data Analysis To analyze the missing data effect on model performance, we studied the correlation between MAPE and missing values. Figure 10 depicts the scatterplots and the Pearson correlation between MAPE and missing values for all sensors in the rollout setting. A slight positive correlation is observed for all the models, from 0.1 for the PSTAINet-ILRB to 0.24 for the PSTAINet-ILB. The points with the highest MAPE for the STAINet, PSTAINet-IB, and PSTAINet- ILRB correspond to the Buriasco and Orbassano sensors, which appear somewhat isolated in the scatterplot. The high MAPEs may be related to the large amount of missing data in recent years and to the pronounced decline observed in their measurements. 15 Salis et al Figure 10: Correlation plots between MAPE and missing values percentage for STAINet, PSTAINet- IB, PSTAINet-ILB, and PSTAINet-ILRB; in all plots, the red line is the bisector and r is Pearson’s correlation coefficient. 7.2.2 Temporal Error Analysis In Figure 11 we reported the MAPE averaged over all sensors for each date of the test set, along with the corresponding standard deviation. PSTAINet-ILB showed the lowest mean MAPE, with slightly higher errors in summer 2022 and summer 2023 – the former affected by a severe drought. STAINet, and PSTAINet-ILRB produced bigger errors in summer and autumn 2022 related to the drop, not well predicted, in groundwater levels in many locations (see Carmagnola in Figure 21h, Cavour in Figure 21k, Fossano in Figures 21m and 21n, and Tarantasca in Figure 21u). This behaviour is also observed, though less accentuated, for the PSTAINet-IB model, which produced slightly larger errors also during spring 2023. Nevertheless, given the limited length of the test set, it is difficult to state seasonal error patterns for the models. It would be necessary to have more observations from additional years as a test set – a condition that will be met in the near future, once data from the 2024 and 2025 field measurement campaigns become available. Figure 11: Average MAPE over all sensors in the rollout setting on the test set. Vertical dashed lines represent the MAPE standard deviation for each week, while vertical dense bold lines represent the corresponding monthly averages and standard deviations. 7.2.3 Trend & Seasonality Analysis We also investigated the association between the trend and seasonality of the groundwater level time series and the error, specifically with MAPE and KGE. In more detail, we decompose each time series with a Seasonal and Trend decomposition using LOESS 6 (STL) [45] assuming each time series Y as a summation of three terms, namely a trend component T , a seasonal component S, and a residual component R, in other words Y = T + S + R. Firstly, the trend component T is estimated using LOESS to capture the overall direction of the series by smoothing short-term fluctuations. Then, on the residual between T and the original Y , the seasonal term S is estimated to extract repeated patterns occurring at the yearly scale. Lastly, the residual component R is determined as the remaining fluctuations not captured by the other two terms 7 . We then adopted metrics to quantify the strength of the seasonal and trend components, specif- ically the Trend Strength S T (Equation 15) and the Seasonal Strength S S (Equation 16) [46]: S T = max 0, 1− V ar(R) V ar(T + R) (15) 6 LOESS stands for LOcally Estimated Scatterplot Smoothing, which is a non-parametric local non-linear regression statistical method. 7 For more technical details, please refer to [45] 16 Salis et al S S = max 0, 1− V ar(R) V ar(S + R) (16) Figures 12 and 13 report the scatterplots and the Pearson correlation between MAPE versus S T , and KGE versus S T , respectively. For all models, a fairly strong positive correlation (ranging from 0.39 to 0.52) is present between MAPE and S T while a negative one with KGE (ranging from −0.52 to −0.25), stating that the stronger the trend of each time series, the lower the performances are. Indeed, a stronger trend component might indicate a nonstationary time series, which is more difficult to model, especially if the nonstationarity is caused by external factors (e.g., drought and climate change). Figure 12: Correlation plots between MAPE and Trend Strength S T for STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB; in all plots, the red line is the bisector and r is Pearson’s correlation coefficient. Figure 13: Correlation plots between KGE and Trend Strength S T for STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB; in all plots, the red line is the bisector and r is Pearson’s correlation coefficient. Conversely, as Figures 14 and 15 show, a positive correlation appeared between KGE and S S (ranging from 0.43 to 0.55), while a low, and partially negative, one with MAPE (ranging from −0.12 to 0.08). This means that models do not show any definite relation between the absolute percentage errors and stronger seasonal patterns. However, stronger seasonalities ease the prediction of the dynamical pattern of the time series (higher KGE) – an expected association given that a stronger seasonality means a clearer repeated pattern through time. Figure 14: Correlation plots between MAPE and Seasonal Strength S S for STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB; in all plots, the red line is the bisector and r is Pearson’s correlation coefficient. 17 Salis et al Figure 15: Correlation plots between KGE and Seasonal Strength S S for STAINet, PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB; in all plots, the red line is the bisector and r is Pearson’s correlation coefficient. 7.2.4 Spatial Error Analysis Figure 16 and Figure 17 report the NBIAS and MAPE for all sensors in the rollout setting, respectively. Overall, the models tended to overestimate groundwater levels, as indicated also by the positive median NBIAS values reported in Table 2. Among all models, the PSTAINet-ILB exhibited the smallest NBIAS, showing a more balanced spatial pattern – with slight underestimations in the northern sensors and overestimations in the southern ones (Figure 16). The Figure 16: NBIAS for all sensors in the rollout setting on the test set. Figure 17: MAPE for all sensors in the rollout setting on the test set. sensor located in La Loggia represented a challenging case for all models, which was underestimated by all of them (as shown also in Figure 21o). This could be due to the high number of missing data in the last year, but also because this sensor is among the few that did not exhibit a pronounced decreasing trend in recent years (see Figure 20 and in particular Figure 20o) – most likely as a result of its proximity to the Po River. Concerning the spatial distribution of the MAPE in Figure 17, a mild tendency can be observed for higher MAPE values in the north-western and south-western sensors across all models (as for the NBIAS) apart from the PSTAINet-ILB, which displays a more even spatial pattern. 18 Salis et al 7.3 Predicted Equation Components Apart from the better performance achieved, the proposed physics-guided models also provided an estimation of the terms of Equation (8). To assess the physical plausibility of the predictions, we analyzed the map predictions of the equation components, in particular, in the rollout setting. Figures 18a, 18b, and 18c show the predicted PDE component on weeks starting on 2022-07-17, 2022-11-23, and 2023-03-12 respectively. PSTAINet-IB predicted a general displacement of ground- water from the mountains to the plain. However, the very high absolute value of ˆ ∆ GW t ∗ and ˆ R t ∗ undermined their physics soundness. Differently, both PSTAINet-ILB and PSTAINet-ILRB provided insightful predictions also from the hydrological point of view. In particular, concerning the diffusion component ˆ ∆ GW t ∗ , PSTAINet-ILB predicted in recharge periods (e.g., 2022-11-13 and 2023-03-12 in Figures 18b, and 18c respectively) a more pronounced groundwater discharge from mountains, while in dryer seasons, given the lower water content and precipitation, the discharge is limited, with water being more retained in the valleys and beneath the mountains. In contrast, the PSTAINet-ILRB model predicted a diffusion dynamic that remained more stable across different periods. Specifically, it captured a general outflow from the valleys toward the plain. The two learning biased models learnt different physics relationships, and even though both appeared sensible, a rigorous evaluation of these processes is difficult, given the scarce knowledge of ground- water processes under mountains. Looking at the sink/source predicted component ˆ R t ∗ , the PSTAINet-ILB model detected a more pronounced and extended recharge areas during recharging seasons (Figures 18b, and 18c), while a more sparse pattern with alternated recharge and loss zones in dry seasons (Figure 18a) – probably because of more abstractions or higher temperatures in summer. Concerning the ˆ R t ∗ predictions made by PSTAINet-ILRB, the recharge zone regularization effect is clearly visible, limiting positive values of ˆ R t ∗ only into the delimited recharge zones (in green in Figures 18a, 18b, and 18c). Although this constraint is hydrologically reasonable, it may have been too restrictive in our case, where the shallow groundwater body is strongly influenced by surface forcing, both human and natural. Finally, regarding the values of the predicted term ˆ D, all models produced plausible values. However, PSTAINet-ILB and PSTAINet-ILRB reported a more sensible spatial pattern with higher values in the valley that gradually attenuated through the plain. In conclusion, the predicted equation component analysis suggests that the absence of supervision on the equation components may be problematic. Introducing regularization terms on ˆ ∆ GW t ∗ and ˆ R t ∗ enabled the models to generate results that are more in line with physics principles. 7.4 Time Series Reconstruction We further evaluated the general ability of PSTAINet-ILB in reproducing the full available time series. Using PSTAINet-ILB, we predicted all the time series within a rollout setting with a forecast horizon of 26 weeks (i.e., 6 months). This is in line with the real situation in our case study, in which the ARPA releases groundwater level data on a semester basis. Table 3 reports the median evaluation metrics of the PSTAINet-ILB over all sensors, and Figure 19 depicts the predicted time series for Cuneo, Racconigi, Scalenghe, and Vottignasco. In Section C, we reported the metrics (Table 9) and the predicted time series (Figure 22) for all sensors separately. The model demonstrated strong capability in reconstructing the full time series, even during intervals with missing data, accurately reproducing seasonal dynamics as well as short-term fluctuations. Table 3: PSTAINet-ILB time series reconstruction median Metrics, in the rollout setting with a 26 time step forecast horizon. NBIAS RMSE [m] MAPE [%]NSEKGE PSTAINet-ILB-0.02780.44900.14040.63140.7547 19 Salis et al (a) (b) 20 Salis et al (c) Figure 18: Predicted PDE components by PSTAINet-IB, PSTAINet-ILB, and PSTAINet-ILRB for the weeks starting on 2022-07-17 (a), 2022-11-13 (b) and 2023-03-12 (c). The first column reports ˆ ∆ GW t ∗ , the second ˆ R t ∗ , and the last ˆ D. Green lines identify the recharge zones. For ˆ ∆ GW t ∗ , a blue pixel means inflow while a red one means outflow; for ˆ R t ∗ , a blue pixel means recharge while red represents a loss. (a) (b) 21 Salis et al (c) (d) Figure 19: PSTAINet-ILB time series reconstruction in the rollout setting with a 26 time step forecast horizon for a) Cuneo, b) Racconigi, c) Scalenghe, and d) Vottignasco. 8 Conclusion In this work, we aimed to develop deep learning models capable of generating predictions at any desired location within the ROI by leveraging spatially sparse in situ measurements (piezometers) together with spatially distributed meteorological data (weather video). Moreover, we investigated different bias strategies, namely inductive and learning biases, to incorporate physical knowledge derived from the groundwater flow equation into the models. To this end, we first developed a purely deep learning model, STAINet, based on the Multi-Head Attention mechanism. We then proposed a restructured variant, PSTAINet-IB, using inductive bias by explicitly predicting the autoregressive (h P,t ∗ −1 ), diffusion (∆ GW t ∗ ), and sink/source (R t ∗ ) com- ponents of the governing equation. Furthermore, we derived the PSTAINet-ILB model by introducing regularization terms into the loss function to explicitly supervise these components (learning bias), and the PSTAINet-ILRB model by further integrating prior information on recharge zones. The PSTAINet-ILB model, trained adopting both inductive and learning bias strategies, achieved the best overall performance, accurately fitting the data while producing physically consistent equation components, thus enhancing model trustworthiness. Although PSTAINet-ILRB also performed well, it achieved slightly lower metrics, probably because of the too restrictive and uncertain nature of the recharge zone prior. The predicted diffusion and sink/source components provide meaningful insights for assessing the physical soundness of the model and offer valuable information on groundwater dynamics in regions lacking piezometer coverage. In particular, while PSTAINet-IB produced implausible ˆ ∆ GW t ∗ and ˆ R t ∗ predictions, both PSTAINet-ILB and PSTAINet-ILRB produced physically consistent predictions, with PSTAINet-ILB especially capturing dynamics evolving accordingly to the recharge seasons. The PSTAINet-ILB model demonstrated similar performance whether leveraging true data as lag inputs or iterating its own predictions over the entire test set (years 2022–2023), thus enabling long- term forecasts all over the ROI. This was further confirmed in the time series reconstruction task, where the complete series for all sensors were predicted using a six-month forecast horizon. In this task as well, the PSTAINet-ILB model achieved remarkable results – particularly given the challenges posed by missing data. In conclusion, the adoption of inductive and learning biases has led to improved performance and enhanced the interpretability of deep learning models, despite several limitations remaining. With respect to the inductive bias, a limitation is the requirement of a closed-form equation that has to be explicitly hard-coded into the model – a condition that cannot always be satisfied in practice, especially in cases of partial or incomplete knowledge. Differently, under the learning-bias strategy, 22 Salis et al the governing equation constraint must be formulated as a loss term. This can substantially increase the computational cost (particularly due to gradient evaluation) when dealing with deep and complex models. Furthermore, the number of control points, namely the locations at which the learning-bias loss terms are evaluated, has a significant impact on both computational cost and the effectiveness of prior injection. On the one hand, increasing the number of control points strengthens adherence to the governing equation; on the other hand, it raises computational demands. This trade-off must therefore be carefully considered and evaluated on a case-by-case basis. Ultimately, even if the weight of the learning bias loss terms could be tuned, a strict compliance of the model’s output with the governing equation cannot be guaranteed. For this reason, adopting both inductive and learning bias strategies jointly may be beneficial to inject prior domain knowledge more rigorously and robustly. The proposed models define a general pipeline applicable to a wide range of environmental appli- cations, where in situ measurements and spatially distributed data are jointly leveraged to predict a target variable, and where a governing equation, even if with some unknown parameters, is desirable to be incorporated. Physics-guided machine learning thus enables the development of hybrid hydro- logical and, more generally, Earth system models that combine the strengths of theory-driven and data-driven approaches, improving computational efficiency while retaining the flexibility to adapt to observed data. In this way, such models provide more effective tools for representing and predicting complex Earth system dynamics. 9 Future Works Even though we obtained remarkable results, several research directions can be undertaken. First, given the widespread presence of missing and irregular data across many disciplines, a more exten- sive evaluation of alternative strategies for handling missing data would be valuable. For instance, attention-based mechanisms could be exploited to design masking strategies that avoid explicit im- putation. In addition, alternative neural architectures could be explored to capture better complex and context-dependent relationships, such as Mixture-of-Experts [47]. In developing the proposed physics-guided strategies, a simple yet effective Euler integration scheme was adopted; however, future work could investigate Neural Differential Equation frame- works to assess more advanced numerical integration methods [48, 25]. While such approaches would increase computational demands, they could enable models to better learn the underlying dy- namics and provide a more rigorous mathematical foundation for incorporating additional physical constraints (e.g., mass conservation). Finally, future studies could benefit from the inclusion of additional data sources, including newly acquired groundwater level observations from the 2024 and 2025 field campaigns, as well as other spatially distributed datasets such as land-cover maps [49] or soil moisture products [50]. These data could offer more informative exogenous forcings, particularly for those related to anthropogenic activities (e.g., irrigation). Author Contributions Matteo Salis: Writing - original draft, review & editing, Investigation, Conceptualization, Visu- alization, Validation, Software, Methodology, Data curation, Resources. Gabriele Sartor: Super- vision, Software, Investigation, Conceptualization, Resources. Rosa Meo: Supervision, Writing - review. Stefano Ferraris: Supervision, Conceptualization, , Writing - review. Abdourrahmane M. Atto: Supervision, Investigation, Conceptualization, Writing - review. Acknowledgments We acknowledge ISCRA for awarding our Deep Learning for Spatio Temporal Phenomena, application to environmental data (DL4STP) project access to the LEONARDO supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CINECA (Italy). The authors declares that they have no conflict of interest to discosle. The authors disclose the use of LLM-based tools for grammar and spelling. Data & Code availability All adopted data are freely accessible, in details the groundwater data at the ARPA website https: //shorturl.at/F3KzQ, and ERA5-land at 10.24381/cds.e2161bac. The code is released in a GitHub repository at https://github.com/Matteo-Salis/physics-guided-gwl with additional animations (GIF) of the predictions over the test set. 23 Salis et al Additional Materials A Sensors Statistics and Time Series In Table 4 we reported summary statistics for every sensor considered in the study. Figure 20 shows the time series plots of the other twenty-four sensors in addition to the four shown in Figure 2. B Test Set Predictions In Table 5, 6, 7, and 8 we reported the evaluation metrics computed on every sensor and every model in the iterative prediction (rollout) scenario. Furthermore, Figures 21 shows models rollout predictions in the test set for the other twenty-four sensors in addition to the four shown in Figure 7. C Time Series Reconstruction In Table 9, we reported the evaluation metrics for the reconstruction task of the PSTAINet-ILB model described in Section 7.4. Figure 22 depicts the weekly prediction from 01-01-2001 for the other twenty-four sensors in addition to the four shown in Figure 19. Table 4: Summary statistics and percentage of missing data for each sensor. Sensor (Municipality - ID)Mean [m] Std. Dev. [m] Missing (%) Bricherasio - 00103510001348.2230.74611.667 Buriasco - 00104110001274.1831.79537.750 Candiolo - 00105110001232.8120.51339.417 Carmagnola - 00105910001231.8930.8133.500 Carmagnola - 00105910002223.7680.4038.333 Cavour - 00107010001274.7771.84222.000 Collegno - 00109010001265.2231.11420.583 La Loggia - 00112710001219.3580.44123.500 Orbassano - 00117110001257.0171.78830.833 Poirino - 00119710001237.1570.54476.917 Scalenghe - 00126010001256.4191.0614.333 Torino - 00127210001227.7300.5644.750 Torino - 00127210003208.7700.35936.083 Virle - Piemonte 00131010001241.8770.58532.500 Barge - 00401210001321.5790.4880.583 Bra - 00402910001275.2020.34122.583 Busca - 00403410001434.6761.86224.917 Caramagna Piemonte - 00404110001251.4980.48853.583 Cavallermaggiore - 00405910001282.7740.50131.167 Cuneo - 00407810001476.6631.8569.833 Fossano - 00408910001351.2090.48514.417 Fossano - 00408910002399.3671.21221.000 Moretta - 00414310002249.2110.34813.083 Racconigi - 00417910001252.4980.7445.833 Savigliano - 00421510001312.8660.26621.833 Scarnafigi - 00421710001282.3211.16931.833 Tarantasca - 00422510001427.3681.4415.917 Vottignasco - 00425010001380.3870.83426.750 24 Salis et al (a) (b) (c) (d) (e) (f) (g) (h) 25 Salis et al (i) (j) (k) (l) (m) (n) (o) (p) 26 Salis et al (q) (r) (s) (t) (u) (v) (w) (x) Figure 20: Groundwater level time series. 27 Salis et al Table 5: STAINet metrics on the test set per sensor in the rollout setting. Municipality SensorNBIAS RMSE[m] MAPE[%]NSEKGE Bricherasio 001035100010.4791.3890.3750.354-0.016 Buriasco 001041100010.3343.1321.0970.1990.431 Candiolo 001051100010.1490.5640.2050.4400.328 Carmagnola 001059100010.2361.0290.4000.4060.293 Carmagnola 001059100020.0850.3300.1220.6250.366 Cavour 001070100010.1712.1750.6570.2100.258 Collegno 001090100010.2831.4720.5150.2680.224 La Loggia 00112710001-0.2140.6240.234-3.188-1.148 Orbassano 001171100010.2303.4711.3110.067-0.801 Poirino 00119710001----- Scalenghe 001260100010.3151.7080.6250.2880.577 Torino 001272100010.0590.5010.182-0.640-0.753 Torino 001272100030.0890.4500.1620.181-0.257 Virle Piemonte 001310100010.0840.6220.1890.5620.435 Barge 00401210001-0.0150.4270.1040.4010.592 Bra 00402910001-0.0110.3600.0980.136-0.116 Busca 00403410001----- Caramagna Piemonte 004041100010.2480.8350.2870.374-0.534 Cavallermaggiore 00405910001----- Cuneo 004078100010.3142.5250.4980.3300.740 Fossano 004089100010.6081.4910.3990.106-0.188 Fossano 004089100020.3922.0340.4500.1440.499 Moretta 004143100020.0300.2790.0920.6440.442 Racconigi 004179100010.1020.7890.2520.4240.210 Savigliano 004215100010.0230.2250.0570.5230.312 Scarnafigi 004217100010.1411.5610.4580.3230.218 Tarantasca 004225100010.4992.9740.6370.1620.031 Vottignasco 004250100010.2341.3530.2850.058-0.094 Table 6: PSTAINet-IB metrics on the test set per sensor in the rollout setting. Municipality SensorNBIAS RMSE[m] MAPE[%]NSEKGE Bricherasio 001035100010.5171.4820.4060.265-0.165 Buriasco 001041100010.2002.0780.6750.6470.074 Candiolo 001051100010.0980.3800.1280.7450.618 Carmagnola 001059100010.1830.8290.3180.6140.497 Carmagnola 001059100020.0170.1830.0660.8850.742 Cavour 001070100010.0471.2690.3910.7310.183 Collegno 001090100010.1720.9340.3140.7050.667 La Loggia 00112710001-0.1880.5060.201-1.755-0.575 Orbassano 001171100010.0861.3720.4910.854-0.392 Poirino 00119710001----- Scalenghe 001260100010.2561.4480.5070.4890.564 Torino 001272100010.0140.3970.142-0.031-0.296 Torino 001272100030.2850.6130.261-0.5190.167 Virle Piemonte 001310100010.0770.4410.1510.7790.746 Barge 00401210001-0.0110.4220.1000.4140.520 Bra 004029100010.0430.2520.0720.5760.437 Busca 00403410001----- Caramagna Piemonte 004041100010.2500.7340.2890.5170.798 Cavallermaggiore 00405910001----- Cuneo 004078100010.2021.8090.3400.6560.767 Fossano 004089100010.3890.9430.2560.6420.582 Fossano 004089100020.3241.6580.3660.4320.568 Moretta 004143100020.0560.2700.0890.6660.604 Racconigi 004179100010.0800.6540.2260.6040.361 Savigliano 00421510001-0.0710.1990.0490.6270.643 Scarnafigi 004217100010.0691.2840.3910.5420.149 Tarantasca 004225100010.2611.7150.3440.7210.063 Vottignasco 004250100010.0470.6450.1420.7860.596 28 Salis et al Table 7: PSTAINet-ILB metrics on the test set per sensor in the rollout setting. Municipality SensorNBIAS RMSE[m] MAPE[%]NSEKGE Bricherasio 001035100010.2760.7950.2160.7880.531 Buriasco 001041100010.0891.2300.3630.8770.403 Candiolo 00105110001-0.1220.4180.1430.6920.567 Carmagnola 001059100010.0210.4440.1540.8890.331 Carmagnola 00105910002-0.0650.2530.0950.7790.675 Cavour 00107010001-0.0371.2340.3280.7460.298 Collegno 00109010001-0.0010.3120.0940.9670.855 La Loggia 00112710001-0.3990.9280.413-8.280-0.045 Orbassano 001171100010.0130.5760.2030.974-0.204 Poirino 00119710001----- Scalenghe 001260100010.0700.6000.1820.9120.745 Torino 00127210001-0.3780.8660.355-3.8990.069 Torino 00127210003-0.0650.3660.1410.459-0.005 Virle Piemonte 00131010001-0.0430.4330.1440.7880.690 Barge 00401210001-0.0790.4540.1090.3240.585 Bra 004029100010.0140.2010.0540.7310.580 Busca 00403410001----- Caramagna Piemonte 004041100010.2360.6900.2730.5730.870 Cavallermaggiore 00405910001----- Cuneo 004078100010.0750.9330.1690.9080.835 Fossano 004089100010.2480.6030.1630.8540.611 Fossano 004089100020.2861.5490.3270.5040.561 Moretta 00414310002-0.0050.1990.0650.8190.630 Racconigi 00417910001-0.0010.4470.1600.8150.468 Savigliano 00421510001-0.1180.2260.0580.5190.689 Scarnafigi 004217100010.0121.1550.3370.6290.188 Tarantasca 004225100010.0340.7530.1550.9460.603 Vottignasco 004250100010.0360.5850.1220.8240.635 Table 8: PSTAINet-ILRB metrics on the test set per sensor in the rollout setting. Municipality SensorNBIAS RMSE[m] MAPE[%]NSEKGE Bricherasio 001035100010.4151.1610.3250.5480.295 Buriasco 001041100010.1601.8890.5440.7090.441 Candiolo 001051100010.0690.3570.1340.7750.496 Carmagnola 001059100010.0530.5240.1950.8460.424 Carmagnola 00105910002-0.0830.2970.1080.6950.552 Cavour 00107010001-0.0141.3170.3680.7100.228 Collegno 001090100010.0610.6780.1930.8440.407 La Loggia 00112710001-0.3450.8400.357-6.607-0.585 Orbassano 001171100010.0991.7770.6180.756-0.995 Poirino 00119710001----- Scalenghe 001260100010.2311.3110.4580.5810.671 Torino 00127210001-0.1880.6080.230-1.411-0.599 Torino 001272100030.1570.5110.174-0.057-0.285 Virle Piemonte 001310100010.0500.4810.1650.7380.702 Barge 00401210001-0.0280.3780.0740.5300.422 Bra 004029100010.0250.2430.0690.6060.258 Busca 00403410001----- Caramagna Piemonte 004041100010.1140.5220.1760.756-0.601 Cavallermaggiore 00405910001----- Cuneo 004078100010.2211.8090.3550.6560.815 Fossano 004089100010.4331.0850.2850.527-0.528 Fossano 004089100020.2951.5600.3370.4970.626 Moretta 00414310002-0.0310.2500.0830.7140.609 Racconigi 004179100010.0390.5460.1830.7240.321 Savigliano 00421510001-0.0030.1690.0450.7320.675 Scarnafigi 004217100010.0681.2250.3620.5830.182 Tarantasca 004225100010.2561.6410.3280.7450.316 Vottignasco 004250100010.1160.7600.1560.7020.618 29 Salis et al (a) (b) (c) (d) (e) (f) (g) (h) 30 Salis et al (i) (j) (k) (l) (m) (n) (o) (p) 31 Salis et al (q) (r) (s) (t) (u) (v) (w) (x) Figure 21: Test set predictions in the rollout setting. 32 Salis et al Table 9: PSTAINet-ILB reconstruction metrics per sensor in the rollout setting. Municipality SensorNBIAS RMSE[m] MAPE[%]NSEKGE Bricherasio 00103510001-0.0300.4510.1030.6510.659 Buriasco 00104110001-0.0350.9200.2640.7370.852 Candiolo 00105110001-0.1230.4330.1480.3110.760 Carmagnola 00105910001-0.0380.4710.1650.6770.727 Carmagnola 00105910002-0.0900.2940.1060.4750.785 Cavour 00107010001-0.0261.1310.3130.6240.765 Collegno 00109010001-0.0150.5770.1630.7380.872 La Loggia 00112710001-0.1180.4880.176-0.2150.368 Orbassano 00117110001-0.0331.2130.3090.5400.724 Poirino 00119710001-0.0300.3390.1070.5480.743 Scalenghe 00126010001-0.0180.4720.1490.8070.807 Torino 00127210001-0.1200.5420.1910.0420.629 Torino 00127210003-0.1250.4240.161-0.3600.375 Virle Piemonte 00131010001-0.0580.3390.1120.6710.820 Barge 00401210001-0.0570.3770.0870.4070.709 Bra 00402910001-0.0250.2570.0660.4420.602 Busca 00403410001-0.0241.1540.1960.6210.792 Caramagna Piemonte 00404110001-0.0210.3110.0960.5930.661 Cavallermaggiore 00405910001-0.0250.2620.0690.7320.769 Cuneo 004078100010.0030.9280.1480.7600.872 Fossano 00408910001-0.0800.3550.0850.4840.671 Fossano 004089100020.0040.7400.1350.6380.750 Moretta 00414310002-0.0160.1750.0550.7490.840 Racconigi 00417910001-0.0230.4470.1410.6510.655 Savigliano 00421510001-0.0790.2000.0520.4230.768 Scarnafigi 00421710001-0.0110.6440.1690.7040.709 Tarantasca 004225100010.0010.7740.1400.7250.856 Vottignasco 004250100010.0150.4210.0840.7520.793 33 Salis et al (a) (b) (c) (d) (e) (f) (g) (h) 34 Salis et al (i) (j) (k) (l) (m) (n) (o) (p) 35 Salis et al (q) (r) (s) (t) (u) (v) (w) (x) Figure 22: Test set predictions in the rollout setting. 36 Salis et al References [1] H. Tao, M. M. Hameed, H. A. Marhoon, M. Zounemat-Kermani, S. Heddam, S. Kim, S. O. Sulaiman, M. L. Tan, Z. Sa’adi, A. D. Mehr, M. F. Allawi, S. Abba, J. M. Zain, M. W. Falah, M. Jamei, N. D. Bokde, M. Bayatvarkeshi, M. Al-Mukhtar, S. K. Bhagat, T. Tiyasha, K. M. Khedher, N. Al-Ansari, S. Shahid, Z. M. Yaseen, Groundwater level prediction using machine learning models: A comprehensive review, Neurocomputing 489 (2022) 271–308. doi: 10.1016/j.neucom.2022.03.014. URL https://linkinghub.elsevier.com/retrieve/pii/S092523122200282X [2] M. G. Zanoni, B. Majone, A. Bellin, A catchment-scale model of river water quality by Machine Learning, Science of The Total Environment 838 (2022) 156377. doi:10.1016/j.scitotenv. 2022.156377. [3] M. Andrychowicz, L. Espeholt, D. Li, S. Merchant, A. Merose, F. Zyda, S. Agrawal, N. Kalch- brenner, Deep Learning for Day Forecasts from Sparse Observations, http://arxiv.org/abs/ 2306.06079 (2023). [4] R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, A. Merose, S. Hoyer, G. Holland, O. Vinyals, J. Stott, A. Pritzel, S. Mohamed, P. Battaglia, GraphCast: Learning skillful medium-range global weather forecasting (2023). doi:10.48550/arXiv.2212.12794. URL http://arxiv.org/abs/2212.12794 [5] A. Karpatne, G. Atluri, J. H. Faghmous, M. Steinbach, A. Banerjee, A. Ganguly, S. Shekhar, N. Samatova, V. Kumar, Theory-Guided Data Science: A New Paradigm for Scientific Discovery from Data, IEEE Transactions on Knowledge and Data Engineering 29 (10) (2017) 2318–2331. doi:10.1109/TKDE.2017.2720168. [6] A. Daw, A. Karpatne, W. D. Watkins, J. S. Read, V. Kumar, Physics-Guided Neural Networks (PGNN): An Application in Lake Temperature Modeling, Chapman and Hall/CRC, 2022, p. 353–372. doi:10.1201/9781003143376-15. URL http://dx.doi.org/10.1201/9781003143376-15 [7] J. Willard, X. Jia, S. Xu, M. Steinbach, V. Kumar, Integrating Scientific Knowledge with Machine Learning for Engineering and Environmental Systems, ACM Comput. Surv. 55 (4) (2022) 66:1–66:37. doi:10.1145/3514228. [8] H. V. Gupta, G. S. Nearing, Debates—the future of hydrological sciences: A (common) path forward? Using models and data to learn: A systems theoretic perspective on the future of hydrological science, Water Resources Research 50 (6) (2014) 5351–5359. doi:10.1002/2013 WR015096. [9] M. P. Clark, M. F. P. Bierkens, L. Samaniego, R. A. Woods, R. Uijlenhoet, K. E. Bennett, V. R. N. Pauwels, X. Cai, A. W. Wood, C. D. Peters-Lidard, The evolution of process-based hydrologic models: historical challenges and the collective quest for physical realism, Hydrology and Earth System Sciences 21 (7) (2017) 3427–3440. doi:10.5194/hess-21-3427-2017. URL https://hess.copernicus.org/articles/21/3427/2017/ [10] D. Bolster, K. R. Roche, V. L. Morales, Recent advances in anomalous transport models for predicting contaminants in natural groundwater systems, Current Opinion in Chemical Engi- neering 26 (2019) 72–80. doi:10.1016/j.coche.2019.09.006. URL https://w.sciencedirect.com/science/article/pii/S2211339819300413 [11] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learn- ing framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707. doi:10.1016/j.jcp.2018 .10.045. [12] M. Reichstein, G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, Prabhat, Deep learning and process understanding for data-driven Earth system science, Nature 566 (7743) (2019) 195–204. doi:10.1038/s41586-019-0912-1. URL https://w.nature.com/articles/s41586-019-0912-1 37 Salis et al [13] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440. doi:10.1038/s42254-021-0 0314-5. URL https://w.nature.com/articles/s42254-021-00314-5 [14] Y. Wu, B. Sicard, S. A. Gadsden, Physics-informed machine learning: A comprehensive review on applications in anomaly detection and condition monitoring, Expert Systems with Applica- tions 255 (2024) 124678. doi:10.1016/j.eswa.2024.124678. [15] E. de B ́ezenac, A. Pajot, P. Gallinari, Deep learning for physical processes: Incorporating prior scientific knowledge*, Journal of Statistical Mechanics: Theory and Experiment 2019 (12) (2019) 124009. doi:10.1088/1742-5468/ab3195. [16] J. Darbon, T. Meng, On some neural network architectures that can represent viscosity so- lutions of certain high dimensional Hamilton–Jacobi partial differential equations, Journal of Computational Physics 425 (2021) 109907. doi:10.1016/j.jcp.2020.109907. [17] S. Cuomo, V. S. Di Cola, F. Giampaolo, G. Rozza, M. Raissi, F. Piccialli, Scientific Machine Learning Through Physics–Informed Neural Networks: Where we are and What’s Next, Journal of Scientific Computing 92 (3) (2022) 88. doi:10.1007/s10915-022-01939-z. [18] A. Arzani, J.-X. Wang, R. M. D’Souza, Uncovering near-wall blood flow from sparse data with physics-informed neural networks, Physics of Fluids 33 (7) (Jul. 2021). doi:10.1063/5.0055 600. [19] G. Kissas, Y. Yang, E. Hwuang, W. R. Witschey, J. A. Detre, P. Perdikaris, Machine learning in cardiovascular flows modeling: Predicting arterial blood pressure from non-invasive 4D flow MRI data using physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 358 (2020) 112623. doi:10.1016/j.cma.2019.112623. [20] C. Banerjee, K. Nguyen, C. Fookes, K. George, Physics-Informed Computer Vision: A Review and Perspectives, ACM Computing Surveys 57 (1) (2025) 1–38. doi:10.1145/3689037. [21] B. A. Toms, E. A. Barnes, I. Ebert-Uphoff, Physically Interpretable Neural Networks for the Geosciences: Applications to Earth System Variability, Journal of Advances in Modeling Earth Systems 12 (9) (2020) e2019MS002002. doi:10.1029/2019MS002002. [22] W.-P. Tsai, D. Feng, M. Pan, H. Beck, K. Lawson, Y. Yang, J. Liu, C. Shen, From calibration to parameter learning: Harnessing the scaling effects of big data in geoscientific modeling, Nature Communications 12 (1) (2021) 5988. doi:10.1038/s41467-021-26107-z. [23] C. Shen, A. P. Appling, P. Gentine, T. Bandai, H. Gupta, A. Tartakovsky, M. Baity-Jesi, F. Fenicia, D. Kifer, L. Li, X. Liu, W. Ren, Y. Zheng, C. J. Harman, M. Clark, M. Farthing, D. Feng, P. Kumar, D. Aboelyazeed, F. Rahmani, Y. Song, H. E. Beck, T. Bindas, D. Dwivedi, K. Fang, M. H ̈oge, C. Rackauckas, B. Mohanty, T. Roy, C. Xu, K. Lawson, Differentiable modelling to unify machine learning and physical models for geosciences, Nature Reviews Earth & Environment 4 (8) (2023) 552–567. doi:10.1038/s43017-023-00450-9. [24] P. Das, A. Posch, N. Barber, M. Hicks, K. Duffy, T. Vandal, D. Singh, K. van Werkhoven, A. R. Ganguly, Hybrid physics-AI outperforms numerical weather prediction for extreme precipitation nowcasting, npj Climate and Atmospheric Science 7 (1) (2024) 282. doi:10.1038/s41612-024 -00834-8. [25] J. Song, Z. Song, P. Ren, N. Benjamin Erichson, M. W. Mahoney, X. S. Li, Forecasting high- dimensional spatio-temporal systems from sparse measurements, Machine Learning: Science and Technology 5 (4) (2024) 045067. doi:10.1088/2632-2153/ad9883. [26] S. Kim, J. Nathaniel, Z. Hou, T. Zheng, P. Gentine, Spatiotemporal upscaling of sparse air-sea pCO2 data via physics-informed transfer learning, Scientific Data 11 (1) (2024) 1098. doi: 10.1038/s41597-024-03959-w. [27] S. Jiang, Y. Zheng, D. Solomatine, Improving AI System Awareness of Geoscience Knowledge: Symbiotic Integration of Physical Approaches and Deep Learning, Geophysical Research Letters 47 (13) (2020) e2020GL088229. doi:10.1029/2020GL088229. 38 Salis et al [28] S. Cuomo, M. De Rosa, F. Giampaolo, S. Izzo, V. Schiano Di Cola, Solving groundwater flow equation using physics-informed neural networks, Computers & Mathematics with Applications 145 (2023) 106–123. doi:10.1016/j.camwa.2023.05.036. [29] D. Secci, V. A. Godoy, J. J. G ́omez-Hern ́andez, Physics-Informed Neural Networks for solving transient unconfined groundwater flow, Computers & Geosciences 182 (2024) 105494. doi: 10.1016/j.cageo.2023.105494. [30] J. D. Hughes, C. D. Langevin, E. R. Banta, Documentation for the MODFLOW 6 framework, Tech. Rep. 6-A57, U.S. Geological Survey (2017). doi:10.3133/tm6A57. [31] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, I. Polo- sukhin, Attention is All you Need, in: Advances in Neural Information Processing Systems, Vol. 30, Curran Associates, Inc., 2017. [32] J. Mu ̃noz Sabater, Era5-land hourly data from 1950 to present, copernicus Climate Change Service (C3S) Climate Data Store (CDS) (2019). doi:10.24381/cds.e2161bac. [33] A. Arnab, M. Dehghani, G. Heigold, C. Sun, M. Lucic, C. Schmid, ViViT: A Video Vision Transformer, in: 2021 IEEE/CVF International Conference on Computer Vision (ICCV), IEEE, 2021, p. 6816–6826. doi:10.1109/ICCV48922.2021.00676. URL https://ieeexplore.ieee.org/document/9710415/ [34] D. Ulyanov, A. Vedaldi, V. Lempitsky, Instance Normalization: The Missing Ingredient for Fast Stylization (2017). arXiv:1607.08022, doi:10.48550/arXiv.1607.08022. [35] A. Wunsch, T. Liesch, G. Cinkus, N. Ravbar, Z. Chen, N. Mazzilli, H. Jourde, N. Goldscheider, Karst spring discharge modeling based on deep learning using spatially distributed input data, Hydrology and Earth System Sciences 26 (9) (2022) 2405–2430. doi:10.5194/hess-26-240 5-2022. URL https://hess.copernicus.org/articles/26/2405/2022/ [36] A. Allen, S. Markou, W. Tebbutt, J. Requeima, W. P. Bruinsma, T. R. Andersson, M. Herzog, N. D. Lane, M. Chantry, J. S. Hosking, R. E. Turner, End-to-end data-driven weather prediction, Nature (2025) 1–3doi:10.1038/s41586-025-08897-0. [37] E. Perez, F. Strub, H. De Vries, V. Dumoulin, A. Courville, FiLM: Visual Reasoning with a General Conditioning Layer, Proceedings of the AAAI Conference on Artificial Intelligence 32 (1) (2018). doi:10.1609/aaai.v32i1.11671. [38] E. Vorontsov, C. Trabelsi, S. Kadoury, C. Pal, On orthogonality and learning recurrent networks with long term dependencies, in: Proceedings of the 34th International Conference on Machine Learning, PMLR, 2017, p. 3570–3578. [39] N. Bansal, X. Chen, Z. Wang, Can We Gain More from Orthogonality Regularizations in Train- ing Deep Networks?, in: Advances in Neural Information Processing Systems, Vol. 31, Curran Associates, Inc., 2018. [40] E. Tadmor, A review of numerical methods for nonlinear partial differential equations, Bulletin of the American Mathematical Society 49 (4) (2012) 507–554. doi:10.1090/s0273-0979-201 2-01379-4. [41] I. Loshchilov, F. Hutter, Decoupled weight decay regularization, in: Proceedings of the Inter- national Conference on Learning Representations (ICLR), 2019. [42] M. Salis, A. M. Atto, S. Ferraris, R. Meo, Time distributed deep learning models for purely exogenous forecasting: Application to water table depth predictions using weather image time series, Environmental Modelling & Software 193 (2025) 106568. doi:https://doi.org/10.1 016/j.envsoft.2025.106568. [43] S. Lee, K.-K. Lee, H. Yoon, Using artificial neural network models for groundwater level forecast- ing and assessment of the relative impacts of influencing factors, Hydrogeology Journal 27 (2) (2019) 567–579. doi:10.1007/s10040-018-1866-3. URL https://doi.org/10.1007/s10040-018-1866-3 39 Salis et al [44] A. Wunsch, T. Liesch, S. Broda, Groundwater level forecasting with artificial neural networks: A comparison of long short-term memory (LSTM), convolutional neural networks (CNNs), and non-linear autoregressive networks with exogenous input (NARX), Hydrology and Earth System Sciences 25 (3) (2021) 1671–1687. doi:10.5194/hess-25-1671-2021. URL https://hess.copernicus.org/articles/25/1671/2021/ [45] R. B. Cleveland, W. S. Cleveland, I. Terpenning, STL: A Seasonal-Trend Decomposition Pro- cedure Based on Loess, Journal of Official Statistics (1990) 3–73. [46] X. Wang, K. Smith, R. Hyndman, Characteristic-Based Clustering for Time Series Data, Data Mining and Knowledge Discovery 13 (3) (2006) 335–364. doi:10.1007/s10618-005-0039-x. [47] W. Fedus, B. Zoph, N. Shazeer, Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity, J. Mach. Learn. Res. 23 (1) (2022) 120:5232–120:5270. [48] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, D. K. Duvenaud, Neural Ordinary Differential Equations, in: Advances in Neural Information Processing Systems, Vol. 31, Curran Associates, Inc., 2018. [49] M. Buchhorn, M. Lesiv, N.-E. Tsendbazar, M. Herold, L. Bertels, B. Smets, Copernicus global land service: Land cover 100 m: Collection 3: epoch 2019: Globe, Zenodo (2020). doi: 10.5281/zenodo.3939050. [50] P. ONeill, S. Chan, E. Njoku, T. Jackson, R. Bindlish, J. Chaubell, Smap enhanced l3 radiometer global daily 9 km ease-grid soil moisture, version 4 (2020). doi:10.5067/NJ34TQ2LFE90. URL http://nsidc.org/data/SPL3SMP_E/versions/4 40