A groundwater level spatiotemporal prediction model based on graph convolutional networks with a long short-term memory

CNN

convolutional neural network

DL

deep learning

DT

decision tree

FC

fully connected

GCN

graph convolutional network

GWL

groundwater level

GNN

graphic neural network

GRU

gated recurrent unit

LSTM

long short-term memory

MAE

mean absolute error

ML

machine learning

RF

random forest

RNN

recurrent neural network

RMSE

root mean square error

Groundwater level (GWL) serves as an essential indicator of the groundwater environment and is influenced by topographical, geomorphological, and meteorologic factors (Scibek & Allen 2006; Zang et al. 2022). Achieving accurate GWL prediction relies on spatiotemporal modeling using rational influential factors and efficient methods (Arabameri et al. 2019). Some traditional and widely used hydrogeological investigation methods, e.g., field surveys, geophysical exploration, well drilling, experimental tests, and numerical simulations, are usually complex, expensive, and time-consuming (Goldman & Neubauer 1994; Qi et al. 2017). Although several groundwater numerical simulation software, e.g., MODFLOW and HYDRUS, include GWL modeling capabilities (Xu et al. 2012; Ghandehari et al. 2024; Zhang et al. 2024c), their practical application is limited by the need for extensive hydrogeological investigation and test data during the time-consuming model development and calibration processes, which is their main limitation in actual application. Additionally, geographical information system (GIS) and remote sensing (RS) techniques have been widely used in GWL modeling (Singh & Katpatal 2017; Qu et al. 2023). However, due to the fluctuating nature of GWL and its nonlinear relationship with multiple influencing factors, these experience-based and linear weighted models that rely on GIS and RS are often subjective and inaccurate in GWL prediction.

Data-driven machine learning (ML) methods aim to approximate the relationship between input influencing factors and output GWL through iterative leaning processes, without the need for explicitly defined physical parameters and models to describe the relationship (Rajaee et al. 2019; Uc-Castillo et al. 2023). Recently, shallow ML methods, e.g., support vector machines (SVMs) (Mackay et al. 2014), artificial neural networks (ANNs) (Gholami et al. 2015), decision trees (DTs) (Wu et al. 2022), random forests (RTs) (Akter et al. 2021; Khan et al. 2023), extreme gradient boosting (XGBoost) (Osman et al. 2021), adaptive neuro-fuzzy inference systems (ANFISs) (Milan et al. 2023; Boo et al. 2024b), and extreme learning machines (ELMs) (Poursaeid et al. 2022) have been increasingly used in GWL prediction to overcome the deficiencies of traditional methods. Although shallow ML models can capture nonlinear patterns in GWL sequences, the long-term prediction results often do not align well with actual observations.

GWL spatiotemporal prediction aims to extract the implicit, unknown, and meaningful relationship between GWL and the spatiotemporal sequence of influencing factors, in order to predict future GWLs based on this relationship. Therefore, it is necessary to consider not only the temporal interaction of different time points but also the interaction of different spatial locations in the spatiotemporal sequence of GWL and influencing factors (Evans et al. 2020; Ali et al. 2021). The spatial autocorrelation and temporal dependency also affect and interact with each other. Spatiotemporal sequence prediction methods can be based on regular grid and irregular graph data structures. Regular grid data, usually consisting of a series of images or transformed into an image-based model (Wang et al. 2022), can be processed using deep learning (DL) models, e.g., convolutional neural networks (CNNs) (Zhang et al. 2024a, 2024b), residual neural networks (ResNet) (Zhang et al. 2018), densely connected convolutional networks (DenseNet) (Xia et al. 2023), fully convolutional networks (FCNs) (Zhang et al. 2020, 2023b), and long short-term memory (LSTM) networks (Wu et al. 2021; Wunsch et al. 2021), to extract spatial characteristics, auxiliary information, and temporal patterns. However, the spatial distribution of realistic observation points is often irregular, making it difficult for regular grid-based CNNs to extract spatial characteristics from scattered points. To address this, Gori et al. (2005) proposed the concept of graphic neural networks (GNNs), using neural networks on non-Euclidean graphs. Kipf & Welling (2017) presented a scalable approach for semi-supervised learning on graph-structured data based on an efficient variant of CNNs, called graph convolutional networks (GCNs), which operate directly on graphs. GCN-based spatiotemporal prediction models have emerged as powerful tools in various fields, e.g., traffic flow prediction and social network analysis. Bai et al. (2023) proposed a spatiotemporal GNN based on gated convolution and topological attention (STGNN-GCTA), which models intricate spatiotemporal correlations in traffic flow with high accuracy and efficiency. Kumar et al. (2023) introduced a dynamic graph convolution LSTM (DyGCN–LSTM), which effectively captures complex spatial and temporal relationships among remote sensors, thereby enhancing traffic forecasting. Skarding et al. (2024) demonstrated that integrating heuristic methods with GCNs significantly improves dynamic link social network predictions. Therefore, to simultaneously capture spatial autocorrelation and temporal dependency, a GCN-based GWL spatiotemporal prediction model should be developed using irregular graphs to effectively model the intricate spatiotemporal dynamics of groundwater systems.

Using easily measurable and widely available influencing factors, GCN-based spatiotemporal methods can help improve the accuracy and efficiency of GWL sequence prediction. In this study, a GCN–LSTM model is proposed for GWL spatiotemporal prediction using GWL temporal sequence data from scattered monitoring wells and auxiliary spatiotemporal characteristics, e.g., temperature, precipitation, elevation, slope, slope aspect, and distance to streams. The model aims to reconstruct the nonlinear relationship between future GWL and historical GWL temporal sequences and influencing factors.

The geological strata in the study area primarily consist of Quaternary Holocene alluvium and lacustrine deposits (Q₄), Quaternary Middle Pleistocene alluvium and lacustrine deposits (Q₂), and Mesozoic Upper Cretaceous red sandstone, calcareous mudstone, and sandy conglomerate (K₂), as presented in Supplementary material, Table S1 and Figure S1. The Quaternary deposits are widely distributed near rivers and primarily comprise gravel, sand, peat, clay, etc. The Upper Cretaceous stratum in the study area was formed in a terrestrial sedimentary environment, consisting of sand-mudstone of near and shallow lacustric facies with gypsum interlayers deposited in salt lakes. The distribution of GWL is influenced by meteorologic and hydrological factors, while geological factors determine the rate and trend of GWL changes. The specific water yield of the aquifer, which varies according to the geological unit, plays a significant role in GWL changes when precipitation infiltrates into the underground to recharge the groundwater. A larger specific water yield results in less GWL rise, and vice versa.

In this study, a total of 1,080 GWL original observations were collected from January of 2003 to December of 2017, with six records per month, from 16 wells in the study area observed by the Hunan Institute of Geological Disaster Investigation and Monitoring (http://hndk.hunan.gov.cn, 30 June, 2022), as presented in Supplementary material, Table S2 and Figure S2. Depths from the ground surface to the GWLs of the 16 wells range from 9.69 to 71.21 m.

The GWL data from 2003 to 2017 were divided into two parts, i.e., a training dataset (864 records) covering the period from 2003 to 2014, and a testing dataset (216 records) covering the period from 2015 to 2017. Subsequently, the GWL temporal sequence data were optimally partitioned into multiple periods using a fixed temporal step length L= 3 for all 16 wells. Each sampling interval occupies 5 days, resulting in an actual temporal span of L being 15 days. Each GWL period, of length L, along with its corresponding auxiliary spatiotemporal characteristics, was trained to predict the GWL for the subsequent observing time, as shown in Supplementary material, Figure S3. Thus, the GWL prediction task at a specific observing time was transformed into a regression problem solved by DL using the data from its previous period.

In GWL prediction, DL methods only require GWL temporal sequences and auxiliary spatiotemporal characteristics as input data, without the need for detailed physical or hydrogeological characteristics. Consequently, DL methods are more effective in predicting dynamic changes and the high uncertainty of GWL. Moreover, DL methods can learn the intrinsic regularities of GWL from fewer training data, without requiring interaction information among multiple influencing factors (Rajaee et al. 2019).

For the prediction of GWL temporal sequences for an individual well, a GWL prediction model was constructed by combining a LSTM layer, a fully connected (FC) layer, and a dropout method. This model utilized historical GWLs and temporal meteorological information, i.e., precipitation and temperature. For the spatiotemporal prediction of GWLs across all monitoring wells, a GCN–LSTM model was proposed to simultaneously capture the spatial autocorrelations and temporal dependencies among GWL sequences and spatiotemporal influencing factors. The model incorporated the spatial location and attribute characteristics (distance from steams, elevation, slope, and aspect) of wells as additional inputs to a GCN for extracting the spatial autocorrelation between wells. The prediction results of GCN–LSTM model were compared with those of the LSTM model, which did not include a spatial autocorrelation extraction module, using a set of evaluation metrics.

Recurrent neural networks (RNNs) have a natural advantage in processing data with sequence changes. Sequence data are taken as input and are recursively processed in the transmission direction of the sequence, with all cyclic units linked in a chain (Hopfield 1982). However, standard RNNs suffer from the vanishing and exploding gradient problems when dealing with long-term dependencies (Rumelhart et al. 1986; Werbos 1990). LSTM, as a special type of RNN, overcomes these issues by introducing more complex gate units in the recursive module (Ghasemlounia et al. 2021). LSTM also has modules similar to RNN chains, but recursive modules have a more complex structure of gate units, as shown in Supplementary material, Figure S4. LSTM can effectively capture both short-term and long-term dependencies in the sequence data, enabling it to process long-term sequence data more efficiently (Hochreiter & Schmidhuber 1997).

Since GWL is a highly temporally dependent sequence data, ordinary DL models cannot capture the temporal dependency within the GWL sequence. However, LSTM can effectively capture the dependency between GWL sequence data. Therefore, LSTM was applied to process GWL and meteorological temporal sequence data in the GWL prediction of an individual well.

Multiple studies have demonstrated that stacked LSTM networks have the potential to generate highly accurate GWL predictions (Dey et al. 2021; Ghasemlounia et al. 2021; Boo et al. 2024a). Four models were designed for the GWL temporal sequence prediction of an individual well. These models consisted of two hidden layers, which included a LSTM layer with a FC layer or two LSTM layers, as shown in Supplementary material, Figure S5. Considering the limited data, non-linearity, and complexity of the GWL dynamic change system, and the large number of parameters in GWL DL prediction models, dropout (Srivastava et al. 2014) was used as a regularization method to overcome the overfitting and reduce training time. Dropout temporarily abandons a portion of neurons in the network during each training iteration according to a certain probability. This prevents the model from over-relying on specific neurons and reduces their cooperative adaptation ability.

Similar to CNNs, GCNs are used for feature extraction. However, GCNs operate on graphs and extract spatial and attribute characteristics from graph to make predictions. For a graph G = (V, E), with N nodes, V represents the set of nodes, which are the objects in the graph, and E represents the set of edges, which are the connections among the objects. If each node V_i has D attribute features, then the attribute features can be represented by an N × D matrix X, and the graph can be represented by an N × N adjacency matrix A. These matrices X and A are the inputs to the GCN, as shown in Supplementary material, Figure S6. Currently, GCNs have been extensively applied in geosciences, e.g., slope deformation prediction (Ma et al. 2021), geochemical anomalies recognition (Guan et al. 2022), bedrock mapping (Zhang et al. 2023a), mineral prospectivity mapping (Zuo & Xu 2023), and three-dimensional geological modeling (Hillier et al. 2021).

In GCN, each node obtains its implicit characteristics from itself and all its neighbors. When there are multiple monitoring wells in the study area, GWL prediction can be regarded as a spatiotemporal sequence prediction. The dynamic GWL changes not only depend on time but also vary based on the locations of the wells. GCN can uncover the influence of spatial autocorrelations among wells to predict more accurate GWLs.

The performance of the GWL prediction models was evaluated using three common accuracy metrics (Khan et al. 2023; Boo et al. 2024a), i.e., coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE).

The models were initially trained using the training dataset and subsequently evaluated on the testing dataset. The average performances of the constructed four models for GWL prediction across all 16 wells were quantitatively evaluated using determination coefficient R², RMSE, and MAE. The evaluation results on the testing dataset are presented in Table 1, which shows that the LSTM + FC model performs the best in terms of R², RMSE, and MAE (R² = 0.577, RMSE = 0.511 m, and MAE = 0.309 m), followed by the LSTM + LSTM model. The superior performance of the LSTM + FC model is likely due to its effective utilization the temporal sequence features captured by the LSTM network, integrated with the FC layer's efficient feature combination and transformation, thereby enhancing prediction accuracy. In contrast, the LSTM + Dropout + LSTM model performs the worst across all evaluation metrics. This maybe because, in the case of a small dataset, adding a Dropout layer can lead to overfitting or loss of crucial information, thereby reducing prediction accuracy.

Supplementary material, Figure S7 depicts the differences between the predicted and real GWL values of the four models from January 2015 to December 2017 across all the wells. Among the four models, the network consisting of a LSTM layer and a FC layer predicted GWL values that closely matching the true values, demonstrating good continuity and smooth transition. The LSTM + LSTM model also achieved good prediction results for GWLs. Comparison of these two models reveals that LSTM can capture the time dependence of long time series features of the well, while its fitting ability is insufficient. However, introducing the FC layer solves this problem and better captures the turning points of GWL peak changes, ensuring good continuity. The FC layer has strong nonlinear fitting capabilities and can optimize LSTM networks with poor fitting, thereby improving the model's prediction accuracy.

LSTM-based models have proven effective for GWL prediction. Sheikh Khozani et al. (2022) combined LSTMs with auto-regressive integrated moving average (ARIMA) models, e.g., DTs and RTs, to capture both linear and nonlinear components of GWL fluctuations. Zeng et al. (2022) assessed lag times using the maximum information coefficient (MIC) algorithm and optimized LSTM hyperparameters with the grey wolf optimizer (GWO) to enhance GWL prediction accuracy. However, since GWLs represent typical spatiotemporal sequence data, it is challenging to account for their spatial autocorrelations using LSTM-based models (Wunsch et al. 2022). The GCN–LSTM model was trained using GWL data from 16 monitoring wells in the study area simultaneously and predicted the dynamic change trend of GWL for all wells from January 2015 to December 2017. Compared with the GWL prediction results from the best-performing LSTM-based model, which achieved R² values ranging from 0.044 to 0.981 at individual wells, the proposed GCN–LSTM model, considering both temporal dependency and spatial autocorrelation, achieved R² values ranging from 0.224 to 0.956. Supplementary material, Figure S8 shows that most of the coefficients R² of GCN–LSTM models are larger than those of LSTM-based models, and most of the RMSEs, and MAEs of GCN–LSTM models are less than those of LSTM-based models. When predicting GWL at same well, the GCN–LSTM spatiotemporal model achieves an average goodness-of-fit approximately 11.21% higher than the LSTM-based temporal model for that specific well.

The fluctuation of GWL is a spatiotemporal process, and the GCN–LSTM model, which considers both spatial autocorrelation and temporal dependency, provides good modeling results for the nonlinear and non-stationary GWL spatiotemporal sequence data. This model demonstrates higher accuracy compared to the LSTM-based model, which only considers the temporal dependency. Therefore, the proposed GCN–LSTM model represents a more comprehensive GWL prediction method, accurately and efficiently reflecting the dynamic changes of GWL in multiple monitoring wells.

This study proposes a multi-well GWL prediction model based on GCN–LSTM, which comprehensively considers the spatiotemporal influencing factors of GWL. The prediction accuracy and reliability of the GCN–LSTM model are verified using the GWL data from 16 monitoring wells in Xiangtan City. The GCN–LSTM GWL prediction model incorporates the temporal dependency and spatial autocorrelation by integrating the spatial autocorrelation and attribute self-similarity features of wells into a graph structure as input. Compared to the temporal sequence prediction model based solely on LSTM, the GCN–LSTM model effectively reflects the spatiotemporal characteristics of GWL dynamic changes in multiple wells, resulting in higher prediction accuracy and better detail in capturing the change trend of GWL in the study area.

Several limitations remain in this study that warrant further investigation. Firstly, the absence of a sensitivity analysis for the GWL spatiotemporal sequence prediction model restricts our understanding of how various factors may influence model predictions. Therefore, future research should prioritize conducting sensitivity analysis to enhance the robustness and applicability of the GWL spatiotemporal sequence model. Additionally, this study does not account for the impact of extreme weather conditions on GWL predictions, especially given the increasing frequency and intensity of such events due to climate change. Future efforts could improve the model's robustness by incorporating the GCN–LSTM framework with climate scenario simulations under various extreme weather conditions, thus evaluating its predictive capabilities in these contexts.

The authors would like to express their gratitude to the MapGIS Laboratory Co-Constructed by the National Engineering Research Centre for Geographic Information System of China and Central South University, for providing MapGIS^® software (Wuhan Zondy Cyber-Tech Co. Ltd, Wuhan, China).

This study was supported by grants from the Hunan Provincial Natural Science Foundation (Grant Nos 2023JJ60188, 2023JJ60190, and 2022JJ30708), the Changsha Municipal Natural Science Foundation (Grant No. kq2208054), the Research Project of Hunan Vocational College of Engineering (GC22YB01), and the National Natural Science Foundation of China (Grant No. 42072326).

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.