A ConvLSTM Conjunction Model for Groundwater Level Forecasting in a Karst Aquifer Considering Connectivity Characteristics

: Groundwater is an important water resource, and groundwater level (GWL) forecasting is a useful tool for supporting the sustainable management of water resources. Existing studies have shown that GWLs can be accurately predicted by combining an artiﬁcial neural network model with meteorological and hydrological factors. However, GWL data are typically geographic spatiotemporal series data, and current studies have considered only the spatial distance factor when predicting GWLs. In karst aquifers, the GWL is affected by the developmental degree of the karst, topographic factors, structural features, and other factors; considering only the spatial distance is not enough, and the real spatial connectivity characteristics need to be considered. Thus, in this paper, we proposed a new method for forecasting GWLs in karst aquifers while considering connectivity characteristics using a neural network prediction model. The connectivity of a karst aquifer was analyzed by a multidimensional feature clustering method based on the distance index and hydrogeological characteristics recorded at observation wells, and a convolutional long short-term memory (ConvLSTM) conjunction model was constructed. The proposed approach was validated through GWL simulations and predictions in karst aquifers in Jinan, China, and four experiments were conducted for comparison. The experimental results show that the proposed method provided the most consistent results with the measured observation well data among the analyzed methods. These ﬁndings demonstrate that the proposed method, which considers connectivity characteristics in karst aquifers, has a higher simulation accuracy than other methods. This method is therefore effective and provides a new idea for the real-time prediction of the GWLs of karst aquifers. at different locations should be considered to improve the prediction of GWLs. This study comprehensively analyzed the connectivity of a karst aquifer based on the distances between observation wells and the hydrogeological characteristics of the aquifer recorded at these wells, and the connectivity analysis results were incorporated into the developed ConvLSTM neural network prediction model. On this basis, four artiﬁcial neural network prediction models were designed to simulate and predict the GWL in the Jinan karst water area. The experimental results show that the multivariable ConvLSTM model that considers the connectivity characteristics among observation wells has a higher simulation accuracy than the other analyzed models. The results of this research show that this method is effective and provides a new idea for obtaining real-time GWL predictions in karst water areas.


Introduction
Groundwater serves as a critical source of water for domestic water uses, agricultural irrigation, and industrial uses, and the rational use of water resources can support sustainable development. In contrast to surface water characteristics, the volume and dynamic evolution characteristics of groundwater are difficult to directly obtain; these features must be revealed by hydrogeological surveys and long-term sequence analyses. Among these characteristics, the groundwater level (GWL) is the most important index and can reveal the influence of human activities, meteorological conditions, and other factors on the groundwater environment [1][2][3]. Therefore, GWL prediction research has become a hot topic.
GWL prediction methods can be divided into methods that use deterministic models and those that use stochastic models. Deterministic models can estimate the hydrological process in a specific region based on the governing equation of the groundwater flow, but they require large quantities of long-term hydrometeorological data and many parameters to describe the physical characteristics of the studied aquifer system, such as MODEFLOW, FEFLOW [4,5]. For regional simulations, a groundwater aquifer is often generalized, causing difficulties in reflecting local details and introducing great errors at some stations. In deterministic models, the GWL prediction process is complex, and the acquisition of parameters is expensive and time consuming [6].
Therefore, in cases of insufficient groundwater aquifer information, many studies have focused on stochastic models to forecast GWLs. Stochastic models are data-driven models that are established by using probability statistics theory. These models do not require sufficient information on hydrogeological aquifer parameters in the model calibration process [6][7][8]; thus, they are more suitable for predicting GWLs when no detailed groundwater attribute data are available.
The stochastic models used for GWL forecasting mainly include time series models and neural network models. Time series models, in which the use of the autoregressive integrated moving average (ARIMA) is typical, have been widely used in GWL prediction research [9]. However, ARIMA is limited because the GWL is affected by many factors that cannot be considered by these models. For example, the GWL is closely related to rainfall, and the ARIMA cannot describe this correlation or its delayed effect [7], although the improved ARIMAX model can be used to consider the relationship between exogenous variables and groundwater level [10]. However, time series models are limited in that their application is restricted to linear processes of a system, so their use is questionable for predicting groundwater levels using time series models, especially in complex environments [11,12].
Currently, with the development of deep learning technology, artificial neural network (ANN) models have been gradually applied in time series data research such as GWL research, and some achievements have been obtained [13,14]. Compared with time series models, ANNs can better identify the nonlinear behaviors of GWL time series. ANNs can learn and summarize inherent regularity from a set of provided examples and do not require information on the interactions among multiple factors that affect the GWL [15]. Common neural network models used for predicting GWLs include the back-propagation (BP) neural network, radial basis function (RBF) neural network, and long short-term memory (LSTM) neural network [16][17][18][19].
In recent years, the LSTM neural network has been developed rapidly [20][21][22]. Long short-term memory (LSTM) neural network, a kind of recurrent neural network, can simulate and predict GWLs using only historical data [21][22][23]. Existing studies have shown that LSTM prediction models have higher efficiencies and accuracies than time series models and BP neural networks [24][25][26].
However, since GWL data are typical geographic spatiotemporal series data, LSTM prediction models consider the temporal autocorrelation of GWL data as well as the spatial correlations among different observation wells in the prediction process [27]. Existing studies have considered only the spatial distance factor and have not involved analyses of spatial correlations among sites under the influence of complex factors [27][28][29]. GWL dynamic characteristics are closely related to external factors, such as rainfall, and internal hydrological characteristics, especially in karst aquifers. Due to the differentiation that occurs during the formation of karst aquifers, the spatial heterogeneity of karst water is very obvious. The flow process within a karst aquifer includes flow in the karst network and flow in the matrix characterizing dual-medium flow [30]. The flow features in a karst aquifer are closely related to its spatial connectivity, which is affected by the developmental degree of the karst, the local topography, structural features, and other factors [31]. Observation wells with relatively nearby spatial distances may have weak mutual connectivity because of the spatial heterogeneity of karst development. In contrast, observation wells with relatively far distances may have strong mutual spatial connectivity. The influence of the spatial connectivity of karst aquifers on GWL forecasting should thus be considered.
Therefore, we propose a new convolutional LSTM conjunction model (ConvLSTM) for forecasting the GWL in a karst aquifer considering its connectivity characteristics. In this study, the connectivity characteristics of the studied aquifer were analyzed based on the multidimensional k-means clustering method, and the ConvLSTM prediction model was constructed based on a convolutional neural network and LSTM networks. Thus, the proposed method considers the temporal autocorrelation of GWL data but also the spatial correlation among different observation wells, and the proposed method was validated through its application to GWL predictions in a karst aquifer in Jinan, China.

Study Area
Jinan, the capital city of Shandong Province, is located in a mid-latitude inland zone and belongs to a warm, temperate continental monsoon climate zone in northern China. The annual average precipitation is 641.68 mm , and the annual average evaporation is 1500-1900 mm [32]. Jinan is famous for its springs and is named "The City of Springs"-108 springs are distributed within 2.6 km 2 in the center of Jinan [33].
Jinan is a typical area containing distributed karst water resources. Karst groundwater is an indispensable and important resource for industrial and agricultural production and social development in Jinan. The carbonate karst aquifers in Jinan are mainly located in the Majiagou Group in the Cambrian-Ordovician strata. Its lithology is dominated by thick microcrystalline limestone, dolomite, and argillaceous dolomite [33,34].
The main supply source of the karst aquifer system in the Jinan spring catchment is meteoric precipitation. The groundwater runoff of the karst aquifer in the spring catchment is controlled by the local topographic features and lithological and geological structure features, and the runoff has obvious spatial heterogeneity [34]. Figure 1 shows the general situation of the study area; the observation wells selected in this paper are mainly located in the carbonate karst aquifer between the Mashan Fault and the Ganggou Fault in Jinan spring catchment.

Data
The data used in this study include long-term GWL observation data, rainfall data, geological structure data, and hydrogeological parameter data.

GWL observation data
The GWL time series data collected in the study area represent the period from 2009 to 2012, and these GWL data were sampled 6 times a month with a sampling interval of 5 days. We selected 16 observation wells from which to obtain experimental data; these wells of the karst aquifer system in the Jinan spring catchment are numbered No. 1~No. 16, and the spatial location distribution of these GWL monitoring wells is shown in Figure 1.

Rainfall
GWL change characteristics are closely related to the atmospheric rainfall process; to obtain rainfall data, we selected 12 rainfall observation stations in the Jinan spring catchment study area that recorded data from 2009 to 2012. The daily rainfall data were resampled 6 times a month with a sampling interval of 5 days in the temporal domain, and the rainfall data and the GWL data recorded in observation wells were matched in the temporal and spatial domains by interpolation and by extraction to points in the spatial domain.

Hydrogeological characteristics data
The hydrogeological characteristics data included the groundwater depth, terrain slope, relief amplitude, variance in the GWL, distance from a fault, and water yield of a single well (Table A1). The terrain slope and relief amplitude represent the topographic variation characteristics of the area where the corresponding observation well is located. The GWL variance was used to indicate the degree of variation in the water level. The distance from a fault indicates the degree of karst development. The water yield of a single well indicates the ability of the aquifer to produce groundwater, which is an important index to classify the water yield grade. In the paper, spatial clustering was firstly carried out using the spatial distances between the observation wells, and then the hydrogeological characteristics were taken as the clustering index set, and another cluster analysis was used to study the connectivity of the karst water.

K-Means Clustering of Multidimensional Features
The k-means algorithm is a popular partitional clustering method, based on the idea of using the cluster centers (means) as representatives of each cluster, which is widely used [35]. The major factors that can impact the performance of clustering algorithms are choosing the initial centroids and estimating the number of clusters [36]. In this paper, the silhouette coefficient was used to select the optimal number of clusters and evaluate the clustering performance, which considers both the intra-cluster and inter-cluster distances for cluster validation [10,37,38]. First, the classification data are divided into k groups with a k-means clustering algorithm, and then the average contour value of the current iteration of each index k is calculated within the range of the predefined minimum and maximum cluster numbers. Finally, the index with the largest average contour value is selected as the optimal number of clusters for the classification data.
In this paper, the k-means clustering algorithm was used to first cluster the observation wells based on their mutual spatial distances; then, further clustering was performed based on the hydrogeological characteristics index values recorded at each station. Through this method, the observation wells with strong spatial connectivity were divided into the same category, while observation wells with weak connectivity were divided into different categories.

LSTM and ConvLSTM
LSTM was originally proposed by Hochreiter and Schmidhuber in 1997 [16], it improved the cyclic neural network model by introducing gates and well-defined storage units to solve the gradient disappearance and gradient explosion problems, respectively, in time series with excessively long processing times in the recurrent neural network model [20]. The neuron structure of the LSTM network is shown in Figure 2. Although an LSTM network can effectively extract the temporal characteristics of time series, the network cannot capture the spatial features of the data. To better integrate the karst water connectivity analysis, this research designed a ConvLSTM module to analyze the temporal and spatial characteristics of GWL changes [39]. The ConvLSTM module includes a convolutional neural network and LSTM network; the resulting module has the temporal modeling ability of an LSTM but can also depict local features such as a convolutional neural network (CNN) [40]. The module captures basic spatial features by performing convolution operations in multidimensional data and replaces the matrix multiplication step with the convolution operation of each gate in the LSTM unit to extract the spatiotemporal characteristics of GWL changes, characterizing the main component of the GWL prediction model developed in this study.
The internal structure ( Figure 3) and calculation process of ConvLSTM (Equation (1)) are as follows: The calculation process of ConvLSTM is as follows: where i, f , c and o are the input gate, forget gate, control unit and output gate in the LSTM structure, respectively, σ is the nonlinear activation function, x t is the input at time t, W xi , , W co are the weight matrix parameters (for example, W xi is the weight matrix value from the input to the input gate), represents the Hadamard product, * represents the convolution operation, h t represents the output value at time t, and o t represents the gating information in the output gate. The ConvLSTM structure used in this article is shown in Figure 4. It consists of two convolutional layers and two LSTM layers and introduces the attention mechanism [41].

Description of the GWL Prediction Process
In this paper, a ConvLSTM conjunction model that considers connectivity characteristics was proposed for GWL forecasting in a karst aquifer. First, the connectivity of the karst water was analyzed by the mixed clustering method, and the observation wells were classified by spatial clustering based on their mutual distances and the results of the attribute clustering based on the hydrogeological characteristics of each well.
Then, the ConvLSTM neural network model was constructed to simulate and predict GWLs; this model has the time series modeling ability of an LSTM network and can also describe local spatial characteristics such as a CNN. When predicting the GWL of each observation well, the rainfall data of the observation well, historical GWL data, and the historical GWL data of associated wells with strong connectivity are selected as relevant variables; the model comprehensively considers the GWLs, meteorological factors, and spatial correlation and heterogeneity between observation wells to simulate and predict the GWL of the karst aquifer of Jinan Spring Catchment.

Input Data Consideration
Existing GWL prediction research has shown that the past steps of the GWL time series and precipitation are the most commonly used variables input into artificial intelligence models to predict GWLs [15]. Based on the actual situation and data availability in the study area, this paper used the lasso regression method to extract important variables from the original data based on regression weights [42]. The GWL and rainfall were selected as input variables among the GWL, rainfall, surface runoff, evaporation, temperature. This conclusion is consistent with existing research conclusions.
In this paper, a mixed clustering method was designed to analyze the connectivity of karst water by the k-means clustering algorithm. First, the observation wells were preliminarily clustered based on their spatial distance index values; then, a further cluster analysis was carried out according to the hydrogeological characteristics of the observation wells by k-means clustering. The most suitable number of clusters was determined by the silhouette coefficient, which can evaluate the performance of clustering results.
The final results were obtained by combining the results of these two methods. For any two observation wells, only if they belong to the same category in the clustering method based on spatial distance, and also in the same category in the clustering method based on hydrogeological characteristics, the two observation wells are finally classified into the same category. Otherwise, they are divided into different categories. This method can avoid two kinds of errors. First, it avoids considering only the distances between wells and ignoring whether the observation wells are truly connected. Second, it avoids the phenomenon in which the clustering results are not spatially clustered when considering only their corresponding hydrogeological characteristics. However, this method may cause the problem of classification redundancy, which needs to be improved in the future. Figure 5a shows the results of clustering based on distance index in which wells were grouped into three categories, and the average silhouette coefficient value is the largest. Figure 5b shows the results of clustering based on hydrogeological characteristics, which shows some partial difference, compared with Figure 5a; for instance, No. 11,No. 12,No. 15,and No. 16 in the result of clustering based on distance, would undoubtedly be classified into the same category since the four observation wells are very close in actual distance. However, in the results of clustering based on hydrogeological characteristics, they are quite different in the water yield of a single well and variance in the GWL, so they are divided into two categories. Furthermore, due to the terrain slope and relief amplitude factors, No.10 is divided into different categories; compared with other observation wells, No. 1 is the same.  Therefore, based on the idea of considering the connectivity of karst aquifers, combined with the spatial distance index and hydrogeological characteristics clustering results, using the cross-merging method, all observation wells were finally divided into seven categories, as shown in Figure 5c and Table 1. When predicting each observation well individually, the GWLs of the remaining observation wells belonging to the same category could be taken as an input variable to consider the spatial connectivity between observation wells and reveal the spatial characteristics of the GWL in the karst aquifer. Thus, the input variables ultimately selected in this study were the GWL, past steps of the GWL time series, rainfall, past steps of the rainfall series, and the GWL sequences of wells with the same category.

Data Set Processing
The GWL and related variables of the prediction well are expressed as a one-dimensional vector as where G p and R p are GWL and rainfall value of the target prediction well, and G cat is the average GWL value of the remaining observation wells that belong to the same category as the predicted well. Then, the one-dimensional vectors of multiple time steps are formed into a twodimensional matrix, which is used to represent the input data within a period of time, i.e., it is used as a time window (Equation (3)). Further, using this time window to slide the complete data sequence of the prediction well from the time series direction, multiple time window data of this prediction well can be obtained (Equation (4)).
where X is data in a time window, and m is the size of the time window; G m p , R m p and G m cat are G p , R p and G cat at the moment m in the time window; P k is data of the target prediction well k and X k n is the time window n of the target prediction well k; x n m is the x at the time m in the time window n.
Finally, the data of 16 observation wells are integrated into an input matrix to enter the ConvLSTM model (Equation (5)). The network can learn the GWL common features of the observation wells in the area and realize the real-time and efficient GWL prediction effect for multiple observation wells in the area with a single prediction model.

Experimental Design
To verify the validity of the ConvLSTM conjunction model that considers connectivity characteristics in predicting the GWL, this paper designed the following four experiments. The training stage was from January 2009 to December 2011, and the testing stage was from January through December 2012. In each batch of training, 20% of the training set was divided for data validation. The parameter settings, input variable, and data set of the experimental models are shown in Tables 2 and 3; The experiment contents are as follows:  Figure 5c and Table 1. Similarly, Experiment 4 used the ConvLSTM model to predict groundwater levels. Table 2 describes the network set and data set for predicting the groundwater levels at the 16 wells in this study. The model parameters were given in a large number of cases. The number of hidden nodes (4,5,6,7,8), kernel size of convolution layer (3,5,7), and learning rate (0.001, 0.003, 0.005) comprised a total of 45 combinations, so one optimal combination was determined by a trial-and-error method in the validation stage [6]. Table 3 describes the best combination of parameters for each experimental model.  Table 3. Input variables and model parameters that were selected for four experiments.

Number of Hidden Nodes
Learning

Results
To examine the applicability of the ConvLSTM conjunction model for GWL forecasting in a karst aquifer while considering connectivity characteristics, we used the root mean square error (RMSE) and Nash-Sutcliffe efficiency coefficient (NSE) to evaluate the simulation results. Table 4 lists the RMSE and NSE values of the four GWL simulation and prediction experiments, and Figure 6 shows histograms of the RMSEs. For the SV-LSTM, MV-LSTM, D-MV-ConvLSTM and M-MV-ConvLSTM models, the average RMSE values were 1.38, 0.75, 0.56 and 0.46, respectively. The experimental results showed that the MV-LSTM model, which considered the rainfall factor, performed better than the SV-LSTM model at each well. The D-MV-ConvLSTM, which considered the GWL distance_cat , performed better than the SV-LSTM and MV-LSTM models. Moreover, the M-MV-ConvLSTM model considered the connectivity between observation wells, and its accuracy was further improved.  Furthermore, it can be seen from the NSE values listed in Table 4  With this model, among all wells, 50% had NSE values above 0.9, and 80% had NSE values above 0.85, indicating that the M-MV-ConvLSTM model, which considered both meteorological factors and the connectivity among observation wells, had the best prediction effect and high credibility. Figure 7 shows the prediction results obtained by using the four models for a selected number of wells; for all wells, the results can be seen in Figure A1. The predicted data are the GWL observation data from January through December 2012, which were sampled 6 times a month with a sampling interval of 5 days. It can be seen from the figure that, in general, the four models and the observed values had the same trend. However, the prediction results of the M-MV-ConvLSTM model were more in line with the actual observation results than the prediction results of the other three models. This was because the M-MV-ConvLSTM model considered the GWL data and meteorological data of the target observation well itself and additionally considered the GWL data of the associated observation wells.  Moreover, compared with the single-structure LSTM model, the ConvLSTM model comprehensively analyzed the spatiotemporal characteristics of the data and revealed a certain improvement in the time-lag problem of the prediction model. From the prediction results of each well shown in Figure 7, it can be seen that the forecast lag problem was greatly improved by the M-MV-ConvLSTM model. Therefore, considering the prediction accuracy and improvement of the lag problem, the M-MV-ConvLSTM prediction model proposed in this paper, which considers the connectivity of observation wells, has good applicability for predicting GWLs.

Influence of Water Level Fluctuations on Model Accuracy
As shown in Table 4 (Table 4); we believe that it is located on the edge of the hill (Figure 1), which is relatively isolated and has poor water yield (water yield of single well <500 m 3 /d) and is vulnerable to external factors. The model cannot accurately reflect the dynamic change of GWL at this well, so the accuracy of the No. 2 has not been effectively improved.

Effect of Spatial Connectivity on Model Accuracy
As shown in Figure 5, when clustering based on distance, the four observation wells No. 11,No. 12,No. 15,and No. 16 were undoubtedly clustered into the same category because of the short distances between these wells. However, the clustering method based on hydrogeological characteristics divided wells No. 11 and No. 12 into category 5 and No. 15 and No. 16 into category 4. Figure 9 displays the prediction results obtained for these four observation wells. We found that although the four observation wells were very close in their spatial locations, No. 11,No. 12,No. 15,and No. 16 obviously have different GWL fluctuation characteristics and show typical spatial heterogeneity. Moreover, the water yield of the four wells is significantly different, and they are in two different levels in the water-rich classification of karst aquifers, as shown in Figure 1. From the results presented in Table 4, and the results of Experiment 4, it can be inferred that the RMSE and NSE values of the M-MV-ConvLSTM model are the best for these four observation wells. The results reveal the effectiveness of the conjunction ConvLSTM method proposed in this paper for predicting the GWLs of karst aquifers while considering connectivity characteristics. The mixed clustering method in this paper may cause the problem of classification redundancy, which needs to be improved in the future. Moreover, comprehensive geographic research usually focuses on various geographic elements to explore the relationships and interactions of elements behind geographic phenomena and processes [43]; the results of this article show that meteorological data are very important for predicting GWL, and it will be a very good research method if these data can be coupled with meteorological forecast models and used in the modeling.

Conclusions
In this paper, a ConvLSTM conjunction model and GWL prediction method that considers connectivity characteristics in a karst aquifer were proposed. The method proposed herein mainly verifies the idea that the spatial connectivity of observation wells at different locations should be considered to improve the prediction of GWLs. This study comprehensively analyzed the connectivity of a karst aquifer based on the distances between observation wells and the hydrogeological characteristics of the aquifer recorded at these wells, and the connectivity analysis results were incorporated into the developed ConvLSTM neural network prediction model. On this basis, four artificial neural network prediction models were designed to simulate and predict the GWL in the Jinan karst water area. The experimental results show that the multivariable ConvLSTM model that considers the connectivity characteristics among observation wells has a higher simulation accuracy than the other analyzed models. The results of this research show that this method is effective and provides a new idea for obtaining real-time GWL predictions in karst water areas.

Data Availability Statement:
The meteorological data and other model simulation data sets used in this study are available from the corresponding author upon reasonable request.
Acknowledgments: Acknowledgement for the data support from Yangtze River Delta Science Data Center, National Earth System Science Data Center, National Science & Technology Infrastructure of China.

Conflicts of Interest:
The authors declare that they have no conflicts of interest.