Groundwater Recharge Prediction Using Linear Regression, Multi-Layer Perception Network, and Deep Learning

Abstract

As the largest freshwater storage in the world, groundwater plays an important role in maintaining ecosystems and helping humans adapt to climate change. However, groundwater dynamics, such as groundwater recharge, cannot be measured directly and is influenced by spatially and temporally complex processes, models are therefore required to capture the dynamics and provide scientific advice for decision-making. This paper developed, estimated and compared the performance of linear regression, multi-layer perception (MLP) and Long Short-Term Memory (LSTM) models in predicting groundwater recharge. The experimental dataset consists of time series of annual recharge from the year 1970 to 2012, based on water table fluctuation estimates from 465 bores in the states of South Australia and Victoria, Australia. We identified the factors that influenced groundwater recharge and found that the correlation between rainfall and groundwater recharge was strongest. The linear regression model had the poorest fitting performance, with the root mean squared error (RMSE) being greater than 0.19 when various proportions of training data were considered. The MLP model outperformed the linear regression in the prediction capability, achieving RMSE = 0.11 when 80% of training data was considered. The LSTM model was found to have the best performance, whose root mean squared errors were less than 0.12 when various proportions of training data were applied. The relative importance of influential predictors was evaluated using the above three models.

1. Introduction

Robust groundwater recharge estimates are a primary requirement for effective management of water resources and sustainable use of groundwater [1], which plays an important role in the sustainable development of regional societies and economies [2,3]. Groundwater recharge is one of the most difficult components of the water balance to estimate since it cannot be directly measured [4,5] and it is influenced by spatially and temporally complex processes. Models are usually required to help stakeholders understand groundwater recharge, identify the key processes influencing the rate of groundwater recharge, and to inform pathways for sustainable water resources management [6].

Previous studies on groundwater recharge have outlined numerous methods for estimating recharge, including chemical tracers [7,8,9,10,11], physical methods [12,13,14,15,16,17] and mathematical approaches [18,19]. The most widely used chemical tracer method is the chloride mass-balance method (CMB) [7,8,9,10], because it is conceptually simple and inexpensive to implement. However, the CMB method cannot estimate the negative component (groundwater evapotranspiration) of net recharge [1], therefore cannot be applied to groundwater discharge regions. Physical methods include water balance estimation (WB) [15,16,17] and the water-table fluctuation method (WTF) [4,12,13,14]. The WB method provides an estimate of net recharge—a combination of recharge to and evapotranspiration from the groundwater, while the water fluctuation method provides an estimate of gross recharge [20]. These methods provide different estimates of recharge, as they include or neglect factors such as runoff, evapotranspiration and changes in soil moisture, and represent either instantaneous estimates or historical averages. In mathematical simulation and statistical methods, linear regression is often adopted for estimating groundwater recharge [7,18,19]. Crosbie et al. [7] used regression kriging to predict regional groundwater recharge across the Australian continent. Global regression equations were applied to data-sparse areas while kriging of regression equation residuals was used for data-dense areas. Fu et al. [18] also used multiple linear regression models together with 71 climate variables and 17 non-climate variables to analyze the groundwater recharge in South Australia. Mathematical simulation models can provide new insights into the factors that affect groundwater recharge.

Multiple linear regression (MLR) is a linear approach for modeling the relationship between input parameters and resulting metrics. In recent studies, this method is applied to model and analyze groundwater recharge. Mogaji et al. [21] estimated and predicted groundwater recharge rates based on the relationship between rainfall and geophysical parameters in the southern part of Perak, Malaysia. Figura et al. [22] predicted groundwater temperatures in several aquifers in Switzerland based on the relationship between observed groundwater and regional air temperature. Ebrahimi et al. [23] simulated groundwater level variations in the Qom plain, Iran based on linear regression, neural network and support vector machine.

Artificial neural networks (ANNs) enhance the expressive ability of the system based on the collection of connected nodes called artificial neurons. Shamshirband et al. [24] proposed a multi-wavelet ANN for forecasting of chlorophyll a concentration. More recently, the ANN method has been used to solve groundwater-related problems [25,26,27]. Mohanty et al. [28] applied the ANN model to the weekly prediction of groundwater levels in various bores based on expert knowledge and statistical analysis. The ANN was used to predict the groundwater level in a swamp forest in Singapore based on rainfall and surrounding reservoir levels [29]. Pasandi et al. [30] applied the ANN to estimating water-table depth in Shibkooh, Iran using various ancillary data, such as aquifer bed elevation and aquifer thickness.

Recent studies have shown that deep learning has broad prospects in groundwater recharge. Deep neural networks have been shown to be suitable for groundwater management. Kong-A-Siou et al. [31] proposed a recurrent multilayer perceptron for predicting the water table level using rainfall and pumping discharge data. Jiang et al. [32] applied a super-resolution convolutional neural network for classifying paleovalleys, which are significant in groundwater exploration as productive aquifers are often formed there. In recent years, computational advances in processing speed and data storage mean that numerically intensive analyses are now possible at large scales and at a declining cost. Models based on machine learning [33] and deep learning [34] are widely used in many fields, such as forest cover projection [34], climate forecasting [35], flood and typhoon forecasting [36]. The typical regression methods used in these studies, such as linear regression, neural networks, and deep learning, have potential for improving the performance in predicting groundwater recharge dynamics.

The paper developed, estimated and compared the performance of linear regression, multi-layer perception (MLP) and Long Short-Term Memory (LSTM) models in predicting groundwater recharge based on water table fluctuation estimates from 465 bores from the year 1970 to 2012. The main purpose of this study was to estimate groundwater recharge in an area straddling the South Australian and Victorian border using the three machine learning methods. Machine learning, especially deep learning approaches have the potential to improve the non-linear expression ability of a system, not only improving the performance of the model, but enhancing the stability of the whole model. In Section 2, we describe the research domain and datasets used in this research. In Section 3, three time series models adopted for predicting regional groundwater recharge are presented. In Section 4, we firstly analyze the correlation coefficients of influential predictors for groundwater recharge estimation. Subsequently, we examine the efficiency of the three models on temporal prediction of groundwater recharge. Finally, we measure the relative importance of influential predictors using the three models. The findings are discussed and summarized in Section 5 and Section 6, respectively.

2. Datasets

2.1. Study Area

Our research area is located in the Otway and Murray Basins in south-eastern South Australia, referred to as the South East (Figure 1). The area is characterized by the tertiary confined sands aquifer known as the Dilwyn Formation covering an area of 29,000 km² [4]. It is overlain by the unconfined Gambier/Murray Group limestone aquifer. The area is relatively flat, with the highest altitude in the north-east part of the region and land generally sloping downward south and west towards the coast. The surface in places is undulating owing to the dune/flat systems from Pleistocene marine transgressions [4].

The region has a typical Mediterranean climate, with hot dry summers and cool wet winters [37]. The highest annual precipitation is in the southern part of the region, gradually decreasing in the inland area. Annual precipitation is less than annual potential evapotranspiration in almost all parts of the study area [38]. Leaney et al. [39] found that based on the karstic features of the landscape, any runoff could quickly infiltrate into the groundwater. Crosbie et al. [4] reported that most of the groundwater recharge in the region occurred during winter due to the lower potential evapotranspiration and higher rainfall at that time.

2.2. Groundwater Recharge and Potential Variable Datasets

In this study, time series data of annual groundwater recharge developed by Crosbie et al. [4] were used. The groundwater level data in the experiment were evaluated from monthly or semi-annual measurements and used to estimate groundwater recharge using the water-table fluctuation (WTF) method. This method provides an estimate of groundwater recharge through the analysis of water-level fluctuations in groundwater observation wells [40,41].

The development of the groundwater recharge data used in the paper was fully described by Crosbie et al. [42] and has been already applied in a few studies, such as [4,20]. It includes recharge data from 465 groundwater bores, mostly located in South Australia with the remaining few in Victoria. Figure 1 shows the location of the study area in Australia and the spatial distributions of the 465 bores. The length of the time series on groundwater recharge varies from 3 to 41 years, based on the length of the groundwater observation records. The dataset is suitable for regional long term average study on groundwater recharge, because more than 70% of bores have an observation record of more than 15 years, and these bores are distributed throughout the study area. The atmospheric demand variables and groundwater extraction dataset were also used for predicting groundwater recharge.

The SILO (Scientific Information for Land Owners) Drill data [43] used in this study consists of 0.05° gridded daily data of atmospheric variables across Australia, shown in

Table 1. Potential variables used in this study.

Data Type	Explanation	Variables	Unit
Spatial-temporal	Regional annual rainfall	Rainfall	mm
	Regional April-October rainfall	Rainfall4-10	mm
	Regional May-September rainfall	Rainfall5-9	mm
	Regional annual actual evaporation	ET	mm
	Regional April-October actual evaporation	ET4-10	mm
	Regional May-September actual evaporation	ET5-9	mm
	Regional annual maximum temperature	Maxtem	°C
	Regional April-October maximum temperature	Maxtem4-10	°C
	Regional May-September maximum temperature	Maxtem5-9	°C
	Regional annual minimum temperature	Mintem	°C
	Regional April-October minimum temperature	Mintem4-10	°C
	Regional May-September minimum temperature	Mintem5-9	°C
	Regional annual Morton actual evapotranspiration	AnnMact	mm
	Regional mean wet-spell length	MeWS	day
	Regional max wet-spell length	MxWS	day
	Regional mean dry-spell length	MeDS	day
	Regional max dry-spell length	MxDS	day
	Regional rainfall (≥1.0mm) days annually	RD	day
	Regional rainfall intensity (Rainfall/RD) annually	RI	mm/day
Temporal	Annual regional groundwater extraction	Extraction	mm

. SILO datasets are constructed from observational records provided by the Bureau of Meteorology (BoM), which have been processed to infill missing data with interpolated values using smoothing splining and kriging algorithms. Rainfall, maximum temperature and minimum temperature are directly measured variables, while actual evaporation and Morton actual evapotranspiration were derived from pan evaporation and other measured variables [43]. The data were all sampled daily from an interpolated dataset provided by SILO [43]. All the data can be accessed at https://legacy.longpaddock.qld.gov.au/silo/about.html.

Actual evaporation was derived from Class-A pan evaporation. A monthly (or seasonal) actual evaporation value was calculated by adding daily pan evaporation in the month (or season). Similarly, annual actual evaporation was calculated by adding monthly actual evaporation values over the corresponding year. The Morton actual evapotranspiration (MAET) was calculated by the complementary relationship between areal potential evapotranspiration (APET) and point potential evapotranspiration (PPET) in [44]. In Morton’s model [45], APET was estimated using modified Priestley-Taylor equation [44] based on psychrometric constant, atmospheric pressure, slope of saturation vapour pressure and net radiation at equilibrium temperature. Furthermore, PPET was estimated by solving energy and vapour transfer equations simultaneously [44] based on air temperature, equilibrium temperature, net radiation at air temperature, saturation vapour pressure and actual vapour pressure respectively.

In this study, the data from 1970 to 2012 were considered. Maximum temperature and minimum temperature were aggregated to a monthly average, and rainfall and evaporation were aggregated to monthly totals. The monthly quantities were used to capture the seasonal features of various predictors and their impacts on groundwater recharge. The study area in the paper was the same as that in Fu et al. [18]. Fu et al. analyzed extreme rainfall variables containing 99th percentile of rainfall, 95th percentile of rainfall and maximum daily rainfall. They found that there was no close relationship between the above three finer scale predictors and groundwater recharge. Therefore, we did not select finer scale predictors, except RD (rainfall days greater than 1 mm) and RI (rainfall intensity) for groundwater recharge analysis. Compared with Fu et al. [18], the novel contribution of the paper was introducing the machine learning and deep learning methods for predicting groundwater recharge. In addition, 465 bores in the states of South Australia and Victoria were used in the paper instead of 426 bores only in the state of South Australia in Fu et al. [18].

The potential influential predictors are shown in Table 1. These are yearly time-series from 1970 to 2012 averaged over the study area. How these yearly predictors were calculated is described as follows. Firstly, daily or monthly rainfall and evaporation values were added up over 12 months (such as Rainfall and ET) or several months (such as Rainfall4-10 or ET5-9) for each year (from 1970 to 2012) in each bore. Daily or monthly maximum and minimum temperature values were averaged across 12 months or several months for each year in each bore. The average and maximum wet-spell as well as dry-spell days were calculated in terms of corresponding daily values for each year in each bore. Subsequently, the regional values for these predictors were averaged across all available bores for each year (for early years, data for some bores are not available) in the case study area. The time series data of influential predictors were obtained using the same method as that for the groundwater recharge. That is, if only 50 bores have groundwater recharge data for a specific year, then regional influential predictors for that year were calculated only based on the locations of those 50 bores. This provides spatial consistency between groundwater recharge values and the corresponding influential predictors. According to [18], the seasonal rainfalls during the winter period (May to September) and summer period (April to October) are critical for groundwater recharge prediction in the South East. The seasonal rainfall and evaporation were included as influential predictors in the analysis, as were estimates of Annual Morton actual evapotranspiration, mean and maximum wet/dry spell-length days, rainfall days and rainfall intensity. Fu et al. [18] analysed the impacts of extreme rainfall variables (containing the 99th percentile of rainfall, and the 95th percentile of rainfall and maximum daily rainfall) on recharge. Their results showed that these predictors and groundwater recharge were not very relevant. Therefore, we did not select many extreme rainfall predictors, except RD and RI for groundwater recharge analysis.

The dataset was developed based on measured groundwater extraction [46] from 2009 to 2013. Groundwater extraction values during the period 1970–2008 were estimated based on the drill date for the bore, and the assumption that the rate of average groundwater extraction for each bore was constant.

3. Methods

3.1. Linear Regression

In statistics, linear regression is a regression analysis approach for modeling the relationship between various influential predictors and a resulting variable. That is, a linear model is established to fit a relationship between the components of the explanatory predictor datasets and resulting variables. The method of least squares is commonly used when linear regression is applied.

Given the dataset ${\left\{y_i,x_{i1},\dots ,x_{ip}\right\}}_{i=1}^n$ ($i$ is the time index, starting from 1 to $n$) with p explanatory predictors (x_i1 to x_ip) and groundwater recharge observations (y_i) containing $n$ statistical values (here $n=43$ in our case as our time series data are available from 1970 to 2012, totaling 43 years), the model assumes that the relationship between the vector of explanatory predictors and resulting predicted groundwater recharge is linear. Thus, the model has the following form:

$$h_{\theta }\left(x_i\right)={\theta }_0+{\theta }_1{x}_{i1}+\dots +{\theta }_p{x}_{ip}=x_i\cdot \theta ,i=1,\dots ,n$$

(1)

where h_θ(x_i) is the estimated resulting variable, x_i (x_i1 to x_ip) is the explanatory predictors, θ₀ is intercept, and θ_p is slope coefficients for each explanatory predictor.

In this paper, the explanatory predictors for the linear model are the 20 potential variables listed in Table 1. The result variable $h_{\theta }\left(x_i\right)$ is the spatially averaged values of groundwater recharge. The whole dataset have 43 groups of data from year 1970 to 2012. That is, $i$ is the data subscript. The linear relationship between explanatory predictors and groundwater recharge values are established when the parameter matrix θ of the model is learned through a training process. The predicted values are calculated from the learned parameter matrix. The fitting and error results are obtained from the difference between predicted groundwater recharge using the linear regression model and observed values of groundwater recharge.

3.2. Multi-Layer Perception Network

Multi-layer perception network (MLP) is a typical representative of feedforward artificial neural networks. It consists of three parts: an input layer, a hidden layer and an output layer. All of nodes in the hidden and output layers are neurons using a non-linear activation function. Furthermore, the hidden layer can be composed of multiple layers of neurons. A supervised learning technique is applied for training the MLP network, called backward propagation. The multiple layered structure and non-linear activation function in the MLP network distinguishes it from linear regressions. Non-linear data relationships can be distinguished by the MLP network.

The MLP network is sometimes referred to as traditional neural network, especially when it only has a single hidden layer. In a MLP network, all the neurons in the hidden and output layers use nonlinear activation functions to simulate the action potential of biological neurons. In this paper, a Rectified Linear Unit (ReLU) function is used for the activation function in all of neurons using Equation (2). Compared with logistic sigmoid function and hyperbolic tan function, the ReLU has better performance. It can solve the problems of gradient explosion and gradient disappearance, and maintains the convergence rate in a stable state.

$$ReLU\left(x\right)=\{\begin{array}{c}x\hspace{1em}if x>0\\ 0\hspace{1em}if x\le 0\end{array}$$

(2)

Since the network is fully connected, each node in one layer connects to every node in the next layer with a certain weight. Therefore, the output $c$ of each neuron is:

$$c=\phi \left({\sum }_i{w}_i{a}_i+b\right)$$

(3)

where $a_i$ and $w_i$ are the inputs and weights of the neuron respectively, b is the bias of the current neuron and $\phi$ is the activation function ReLU.

In this paper, the MLP network is composed of three layers (an input, one hidden layer and an output layer) of nonlinearly-activating modes (Figure 2). The input layer has 20 nodes, which corresponds to the 20 potential influential predictors that are listed in Table 1. There are 100 neurons in the hidden layer. All of the weights and biases of the linkages between the potential variables and groundwater recharge values are updated when the parameters of the network are learned based on the training data. Then the predicted values are calculated based on the learned weights and biases when the input data is applied. The output layer only had one neuron, which represented the values of groundwater recharge.

At the beginning of the training, network weights and biases were assigned randomly. The solver for weight optimization is a stochastic gradient-based optimizer (Adam). The algorithm then projected these forward from the input layer to the hidden layer based on Equation (3). The results from the hidden layer were propagated to the output layer and the error between value of the output layer and the observed groundwater recharge in training data was obtained. By propagating the error iteratively back to the network, the connection weights and biases were automatically adjusted until the network error reaches a pre-determined value. The predicted values could then be calculated using the input values of 20 potential variables based on the trained network parameters.

3.3. LSTM Model for Regression

To solve the problem of vanishing gradient in traditional recurrent neural networks, Hochreiter and Schmidhuber [47] proposed LSTM, which introduces a new structure called a memory cell to develop persistent long-term dependencies.

The structure of the whole algorithm is composed of one input layer, one output layer and three LSTM layers. The process of training and predicting depends on the temporal sequence features of the data being predicted. The system has $T$ input vectors, $x^{T=1,\dots ,t_T}$, which represent $t_T$-year groundwater values. Similarly, the system output also has $T$ prediction vectors $x^{T=1,\dots ,t_P}$ indicating $t_P$-year prediction results. For this study, the temporal sequence features from the regional average groundwater recharge was developed based on the 465 bores located in the study area from 1970 to 2012.

In order to avoid the proposed model from being overfitted, a method called dropout was used, with dropout rate 0.25 applied to the output of each LSTM layer. It was a regularization technique for reducing overfitting in the model by preventing complex co-adaptations on training data, and improved the accuracy of prediction data. The mean squared error (Equation (4)) was selected as the loss function to regularize the training process, where ${\widehat{y}}^t$ and $y^t$ represent for the predicted values and observed values in year $t$.

$$loss=\frac{1}{T}{\displaystyle \sum }_{t=1}^T{\left({\widehat{y}}^t-y^t\right)}^2$$

(4)

The internal structure of a memory cell in a LSTM layer is shown in Figure 3. The memory cell is composed a forget gate $f_{l,t}$, an input gate $i_{l,t}$, a new memory unit $n_{l,t}$ and an output gate $o_{l,t}$, where $l$ and $t$ are the number of LSTM layers and the current time step respectively. $W_l^t$ and $U_{l,t}$ are the weight matrices corresponding to the input $x_{l,t}^k$ and previous hidden state $h_{l,t-1}^k$, and $b_{l,t}$ is the bias vectors on the memory cell. The superscripts f, i, n and o represent the forget gate, input gate, new memory unit and output gate respectively.

The forget gate $f_{l,t}$ is calculated as follows:

$$f_{l,t}=\sigma \left(W_{l,t}^f\cdot x_{l,t}^k+U_{l,t}^f\cdot h_{l,t-1}^k+b_{l,t}^f\right)$$

(5)

The logistic sigmoid function $\frac{1}{1+e^{-x}}$ is applied as the activate function, $\sigma$, on the gates.

Similarly, the input gate and output gate are defined as follows:

$$i_{l,t}=\sigma \left(W_{l,t}^i\cdot x_{l,t}^k+U_{l,t}^i\cdot h_{l,t-1}^k+b_{l,t}^i\right)$$

(6)

$$o_{l,t}=\sigma \left(W_{l,t}^o\cdot x_{l,t}^k+U_{l,t}^o\cdot h_{l,t-1}^k+b_{l,t}^o\right)$$

(7)

The hidden state $h_{l,t}^k$ and output cell state $c_{l,t}^k$ are denoted as follows:

$$h_{l,t}^k=tanh(c_{l,t}^k)\ast o_{l,t}$$

(8)

$$c_{l,t}^k=i_{l,t}\ast n_{l,t}+f_{l,t}\ast c_{l,t-1}^k$$

(9)

where $\ast$ is the Hadamard product of two vectors. The new memory unit, $n_{l,t}$, can be calculated as follows:

$$n_{l,t}=tanh(W_{l,t}^n\cdot x_{l,t}^k+U_{l,t}^n\cdot h_{l,t-1}^k+b_{l,t}^n)$$

(10)

The three LSTM layers are followed by an output layer in the fully connected structure. The prediction vector, ${\widehat{y}}^t$, in time, $t$, can be calculated using Equation (11).

$${\widehat{y}}_{}^t=\sigma \left(W^{out}\cdot {\widehat{y}}_{}^{k,t}+b^{out}\right)$$

(11)

where $W^{out}$ and $b^{out}$ are the weight matrix and bias vector in the output layer respectively. A logistic sigmoid function is used as the activate function in the output layer for predicting groundwater recharge in year $t$.

3.4. Model Testing and Comparison

The linear regression, MLP and LSTM models were run three times, trained with 70%, 80% and 90% of the observation data, then their predictions were compared with 30%, 20% and 10% of the observation data as a model validation. The statistical measures used for comparing the model performance were the root mean squared error (RMSE) and coefficient of determination ($R^2$). The RMSE and $R^2$ were calculated between the predicted groundwater recharge values and observed recharge values using Equations (12) and (13), to measure the precision and bias of groundwater recharge predictions, respectively.

$$RMSE=\sqrt{{\displaystyle \sum }_{i=1}^n\frac{{(y_i-{\widehat{y}}_i)}^2}{n}}$$

(12)

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{{{\displaystyle \sum }}_{i=1}^n{(y_i-\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{y}_i)}^2}$$

(13)

where $y_i$ and ${\widehat{y}}_i$ are the observed and estimated values respectively and n is the number of the observed values in the testing data. The $R^2$ should be close to 1 to indicate strong model performance, and the RMSE should be as close to zero as possible. The two measures are applied for estimating the prediction of groundwater recharge by the three models. In the statistical analysis, R-squared is a relative measure of fit. Its advantage is that it can be used to measure the whole trend of fitted data. As the square root of a variance, RMSE is an absolute measure of fit. It can measure the error between the whole prediction data and observation data directly. Therefore, in order to ensure the reliability of the results, the relative indicator $R^2$ and absolute indicator RMSE are used to measure the fitted and prediction effects together.

The linear regression model and the MLP model were implemented based on the Python packages StatsModels and scikit-learn, respectively. The LSTM model was built using Tensorflow, which is an end-to-end machine learning platform in Python. All the modelling, data analysis and visualization were conducted in the environment of Python 3.5. TensorFlow is an open source machine learning library based on data flow graphs. It has better computational graph visualizations, and is supported by Google. As a popular framework, Tensorflow is adopted in the work. The environment Python is selected in the work owing to its cross-platform processing capability, open source code, presence of third-party modules and extensive support libraries.

4. Results

The annual time series of regional groundwater recharge was established by averaging the groundwater recharge estimates temporally for each year and spatially based on data from all available bores. The three models were then developed to predict groundwater recharge using the full set of influential predictors. The relative importance of each predictor in predicting groundwater recharge was measured under the three temporal models. Furthermore, all of weights and biases were initialized using Xavier initializer. Then, all of the parameters were optimized based on a stochastic gradient-based optimizer (Adam). The method of initializing and optimizing parameters can ensure that all of parameters maintain the good performance in training and prediction.

4.1. Correlation Coefficients between Potential Variables and Groundwater Recharge

The time series of annual regional groundwater recharge and potential predictors from 1970 to 2012 were applied to establish the temporal regression models. The Pearson correlation coefficients between regional groundwater recharge and each influential predictor are shown in Figure 4.

As expected, the correlation coefficient between rainfall and groundwater recharge was the highest. Annual rainfall had a correlation coefficient of 0.8, rainfall from May to September had a correlation coefficient of 0.88 and for rainfall from April to October it was 0.9. The rainfall and minimum temperature had a positive correlation with groundwater recharge, while in contrast, the evaporation and maximum temperature had negative correlation with groundwater recharge. This could be because these influential predictors influence the negative component of groundwater recharge, or evapotranspiration. Almost all the absolute values of correlation coefficients for evaporation, maximum temperature and minimum temperature were between 0.3 and 0.4.

Other influential predictors, MeWS, MxWS, RD and RI had positive correlation coefficients, as they are associated with rainfall and infiltration, while MeDS, MxDS and Extraction had negative correlation coefficients as they are associated with evapotranspiration processes. The RD correlation magnitude with groundwater recharge (r = 0.76) was higher than that of RI (r = 0.45). It indicated that groundwater recharge was a process of gradual accumulation, because at an annual time scale, the number of rainfall days was more important than rainfall intensity (Figure 4). As expected, the correlation magnitudes for mean and max Wet-Spell length (r = 0.66, r = 0.63) were higher than that of the mean and max Dry-Spell length (r = −0.50, r = −0.15), as wet days are likely to result in groundwater recharge than dry days. Groundwater extraction had a negative correlation coefficient (r = −0.43) with groundwater recharge. This can be attributed to the fact that groundwater extraction is simply a negative recharge. Furthermore, it has been shown that lower water tables can lead to lower recharge [42], and this is consistent with the negative correlation relationship between groundwater extraction and recharge.

4.2. Temporal Prediction of Groundwater Recharge

The time series of regional groundwater recharge based on the average values of groundwater recharge from 1970 to 2012 are shown in Figure 5a–c, marked as “Observation”. Figure 5a shows the results from the linear regression method, which was trained with the first 80% of the observed data, and then predicted groundwater recharge for the last 20% of the time series. Similarly, Figure 5b,c demonstrate the results from the MLP model and the LSTM model, respectively. Both models were trained using the first 80% of the observed data and predicted the rest 20% of the time series.

As shown in Figure 5a, there was relatively good fitness and several small mismatches between the trained data and observation for the linear regression model. However, the fitting performance between the predicted data and observation was poor, especially from 2010 to 2012, and the mismatch was serious. That is, the time series’ features obtained from the trained data could not successfully be applied to the predicted data when the linear regression model was used.

For the MLP model (Figure 5b), the fitness between trained data and observations were better than that of the linear regression model. The trends in trained data and observations were consistent. However, the fitting result between the predicted data and observations was only slightly better than that of the linear regression model. That is, the MLP method accurately learned the time series characteristics of the trained data, but it still lacked in aspects of data prediction. This can be attributed to the over-fitting of learning from trained data based on the MLP method.

Using the LSTM model (Figure 5c), the fitting performance between trained data and observations was worse than that of MLP model and similar to that of the linear regression model. However, the fitting result between the predicted data and observation was better than those from the linear regression and MLP models, especially from 2010 to 2012. The trends between the predicted data and observation were consistent with those in the LSTM model (Figure 5c), but not consistent with those in the linear regression and the MLP model. Compared with the insufficient learning of the linear regression model and the over-fit learning of the MLP model, the time series’ features learned from the trained data by the LSTM model could be successfully applied for the prediction of groundwater recharge.

The linear regression, MLP and LSTM methods were then compared for different percentages of the observation data being used for training and testing the models (

Table 2. The trained and predicted groundwater recharge results based on the three models and the percentage of observation data used to train the models.

Training Data (%)	Linear Regression			MLP Model			LSTM Model
	$RMS{E}_{Tra}$	$RMS{E}_{Pre}$	$RMS{E}_{Sum}$	$RMS{E}_{Tra}$	$RMS{E}_{Pre}$	$RMS{E}_{Sum}$	$RMS{E}_{Tra}$	$RMS{E}_{Pre}$	$RMS{E}_{Sum}$
70%	0.06	0.20	0.13	0.02	0.19	0.12	0.06	0.12	0.09
80%	0.06	0.19	0.11	0.02	0.11	0.06	0.07	0.10	0.08
90%	0.06	0.20	0.09	0.02	0.13	0.05	0.04	0.11	0.06
Training Data (%)	Linear Regression			MLP Model			LSTM Model
	$R^2{}_{Tra}$	$R^2{}_{Pre}$	$R^2{}_{Sum}$	$R^2{}_{Tra}$	$R^2{}_{Pre}$	$R^2{}_{Sum}$	$R^2{}_{Tra}$	$R^2{}_{Pre}$	$R^2{}_{Sum}$
70%	0.96	0.46	0.79	0.99	0.49	0.82	0.94	0.77	0.88
80%	0.95	0.44	0.85	0.99	0.82	0.95	0.93	0.84	0.92
90%	0.95	−0.11	0.89	0.99	0.54	0.96	0.98	0.70	0.96

). In order to compare the three methods, an error indicator ($RMSE$) and a fitness indicator ($R^2$) were applied into testing the whole prediction results. $RMS{E}_{Tra}$ and $R^2{}_{Tra}$ are the error and fitness values for the training data, $RMS{E}_{\mathrm{Pr}e}$ and $R^2{}_{\mathrm{Pr}e}$ are the error and fitness values for the prediction data, and $RMS{E}_{Sum}$ and $R^2{}_{Sum}$ are the error and fitness values for the total data containing both the training data and the prediction data. All of the data for the influential predictors and groundwater recharge observations were divided into two groups: training data and prediction data. In the experiment, the ratio of training data was selected as 70%, 80% and 90% respectively.

For all the three models, $RMSE$ for all the total dataset ranged from 0.02 to 0.20. For each model, the $RMSE$ for all the data decreased and the fitness indicator $R^2$ increased when the ratio of training data to prediction data increased. It also demonstrates that the prediction error for the training dataset was lower than that for the prediction dataset. The errors for the total dataset produced from the MLP and the LSTM were similar, and the performance of the linear regression model was slightly inferior to the two more complex methods.

The three models all presented excellent fitting performance on the training dataset. All fitting errors (represented by $RMSE$) were lower than 0.07 and all fitness values (represented by $R^2$) were greater than 0.9, no matter which ratio of training data was chosen. For the three models, the fitting performance on the training dataset was obviously better than that on the prediction dataset. Furthermore, compared with the other two models, LSTM had the minimum differences in both fitting errors and fitness values between the prediction dataset and the training dataset. This demonstrated that LSTM model has the higher generalization performance.

Fitting error of the predicted data was the most important factor for comparing the performance of various regression models. For the linear regression, $RMSE$ was greater than 0.19 when the 70%, 80% and 90% ratios of training data were considered, and the fitness $R^2$ was negative when 90% of the data was used to train the model. The negative value represented that the fitting performance was worse than just fitting a horizontal line and linear relationship was not suitable for fitting the data.

For the MLP model, $RMSE$ and $R^2$ were better than that those from the linear regression, with the fitting error lower than 0.20 for the 70% ratio of training data. The fitting performance was satisfactory with fitting error $RMSE$ equalling 0.11 when the 80% ratio of training data was considered. The $RMSE$ values were lower than those of the other two methods for all the three ratios of training data. Also, all the fitness values were more than 0.7, no matter which ratio of training data was chosen. That is, the stability of the LSTM model was the best in the three compared methods.

The LSTM model had the best prediction performance in terms of its $RMSE$ values and $R^2$ values regardless of the proportion of training data chosen. Although the performance of the MLP model was worse than that of the LSTM model, it still gave a better estimate of groundwater recharge than the linear regression method for all of three ratios of training data. It demonstrated that compared with the traditional regression method, deep learning and machine learning methods can greatly improve the whole performance of groundwater recharge prediction.

4.3. Relative Importance of Influential Predictors

The relative importance of potential influential predictors for predicting groundwater recharge was investigated for each of the linear regression, MLP and LSTM models. The coefficients of determination and RMSEs are shown in Figure 6, where influential predictors were excluded one by one from linear regression, MLP and LSTM models. The variable importance ranking was found to be the same whether it was determined through RMSE (Figure 6a) or coefficients of determination (Figure 6b).

For the linear regression model, the least important variables were found to be the minimum temperature (Mintem, Mintem4-10 and Mintem5-9) and Average annual Morton actual evapotranspiration (AnnMact), and the most important predictors were mean dry spell and wet spell lengths (MeDS and MeWS) and annual number of rainfall days (RD), in that order. Rainfall data (Rainfall, Rainfall4-10 and Rainfall5-9) also had an obvious influence. This was unsurprising, because correlations between groundwater recharge and variables of minimum temperature and annual Morton actual evapotranspiration were relatively low (See Figure 4), and the correlations in mean wet/dry spell-length days and rainfall days were significantly high (See Figure 4). Furthermore, it can be attributed to linear regression model which applies linear relationship between potential predictors and is consistent with correlation coefficients between corresponding predictors and groundwater recharge.

For the MLP model, rainfall, maximum wet spell length (MxWS), maximum temperatures (Maxtem, Maxtem5-9, Maxtem4-10), and evaporation (ET, ET4-10 and ET5-9) were relatively insignificant variables, while Mintem, Mintem4-10, Mintem5-9 and AnnMact were found the most important predictors.

For the LSTM model, the least important variables were Maxtem, Maxtem5-9, Maxtem4-10 and RD. For this model, the most important predictors were AnnMact, ET4-10, Mintem, Mintem5-9, MeDS, MeWS, Rainfall4-10 and RI. Overall, the RMSEs were lower for the LSTM model, indicating that it retained its good performance even when a variable had been removed from the inputs.

The relative importance of influential predictors are different for the MLP and LSTM models. The phenomenon can be attributed to two causes. Firstly, over-fitting can be a problem, especially for the MLP model. This model relies entirely on the training dataset to develop its prediction algorithms, which results in its poor performance for the predicted dataset. The generalization ability of the whole model was therefore reduced and the accuracy of the predictor importance measurement was also affected. Secondly, the problem of non-convex optimization makes convergence to a local optimum possible when the parameters of the model are learned and adjusted. This will lead to the deviation of the predictor importance’s estimation. Finally, different local optimal solutions can lead to differences in the relative importance of same influential predictor when MLP and LSTM models are applied respectively.

5. Discussion

In this paper, linear regression, MLP and LSTM models have been used to predict a time series of annual average groundwater recharge from 1970 to 2012, based on 465 bores in the South East of South Australia and Victoria. The performance on the three models in predicting groundwater recharge was assessed using various ratios of trained dataset. The study obtained the following findings.

5.1. Performance and Comparison of Models

Compared with the MLP and linear regression models, the LSTM model showed the best performance in predicting groundwater recharge in the case study area.

Previous research on using machine learning to understand groundwater processes includes the use of a MLP model to predict weekly groundwater levels in the Mahanadi Delta, India [28] and a linear regression model to predict groundwater recharge based on 71 climate variables and 17 non-climate variables also in the South East of South Australia and Victoria [18]. The purpose of this paper is to find the applicable model for groundwater recharge prediction by comparing the linear regression model, MLP and LSTM models based on the same data format. According to Fu’s research [18] in the same study area, seasonal rainfall from May to September and from April to October are the most important influential predictors for groundwater recharge based on multivariate linear regression model. Therefore, we calculate the monthly average quantities based on daily data on various influential predictors from 1970 to 2012, such as maximum temperature, minimum temperature, rainfall and evaporation. The LSTM model has been used to predict groundwater heads [48], but until this paper, it has not been used to predict groundwater recharge timeseries. Results from this study (Section 4) show that the LSTM model consistently obtained the best performance for annual groundwater recharge prediction compared with MLP and linear regression models.

Research in the field of thermodynamics, energy and fuels also shows that the LSTM model performs better than MLP network in forecasting aggregated power load and photovoltaic (PV) power output [49].

The LSTM model is a more complex type of machine learning model that includes three control units: input gates, output gates, and forget gates. As the information enters the model, the control units in the LSTM assess each piece of information for which will be retained and which will be forgotten. The LSTM model can therefore acquire long-term dependencies to solve the problem of gradient vanish and improve predictions of groundwater recharge.

The performance of the MLP model was also found to be better than the more simple linear regression model approach. The MLP model is also known as a feedforward neural network which has one or more hidden layers between the input layer and the output layer of the network. The trained parameters in these hidden layers improve the predictive ability of the model.

An important finding was that the performance of the LSTM model appeared to be almost unaffected for various ratios of trained and prediction dataset were considered. A benefit of this model is that it may reduce the length of time series data required to train the model, and that its results were relatively resilient to changes in parameter settings.

In addition, we found that the model fitness declined when predictors that have low correlations with groundwater recharge were added into the linear regression model. The fitness of the linear regression model depended on the correlation strength between selected predictors and groundwater recharge. That is, the model had better fitness when those predictors with higher correlation coefficients were added into the model. This reflects that the linear regression model was not able to deal with the non-linear relationship among the predictors. On the contrary, different combinations of influential predictors had little influence on the fitness of the MLP and the LSTM. That is, compared with the linear regression model, the two models had better robustness in predicting complex data relationship. Therefore, the hidden layer in MLP model and the memory unit in LSTM model captured the non-linear relationship among influential predictors, and improved the robustness of the models.

Traditional time series analysis methods, such as the autoregressive integrated moving average with explanatory variable (ARIMAX) model, are difficult to capture non-linear relationship. The hidden layer in MLP model and the memory unit in LSTM model can help improve the non-linear learning ability of the models.

5.2. Influential Predictors Identification

The relative importance of all predictors to drive the temporal changes in groundwater recharge was measured using the linear regression, MLP and LSTM models. The results showed that the relative importance of influential predictors was different for each of the three models.

The most important predictors for the linear regression model were the mean number of days wet and dry, number of rain days and three rainfall metrics (MeDS, MeWS, RD, Rainfall, Rainfall4-10 and Rainfall5-9). This is consistent with standard hydrogeological understanding of groundwater recharge, which often estimates recharge as a proportion of rainfall or difference between rainfall and evapotranspiration [22].

For the MLP model, the important predictors were almost opposite to the hydrogeological understanding of recharge: minimum temperature and annual evaporation (Mintem, Mintem4-10, Mintem5-9 and AnnMact), although evaporation is known to decrease net groundwater recharge [42]. Future research is advisable to determine whether these relationships were coincidental or are consistent across many different locations and climate types.

For the LSTM model, the most important predictors were found to be a combination of evapotranspiration, minimum temperature, annual and seasonal evaporation, mean number of dry and wet days, rainfall and rainfall intensity (AnnMact, ET4-10, Mintem, Mintem5-9, MeDS, MeWS, Rainfall4-10 and RI). The LSTM model had a consistently high performance when any one of the influential predictors was removed from the analysis and it was less sensitive to these changes in input. This model appeared to rely on more of the influential predictors to predict groundwater recharge, therefore was more robust when any of the predictors was removed from the analysis. Again, future research is required to determine how consistent these important influential predictors are across multiple locations and climates.

5.3. Implications for Groundwater Management

Accurate evaluation of groundwater recharge is the foundation for equitable and reliable water resources planning. This research verified that machine learning and deep learning methods can be applied to the estimation of groundwater recharge. Out of the three models, the LSTM model had the most precise groundwater recharge predictions. The results from this study provide groundwater researchers and managers a valid reference for the selection of appropriate machine learning models in the future.

This research has explored the groundwater recharge estimation and important drivers for changes in temporal groundwater recharge in the South East of South Australia and Victoria. The performance of the models in this study may also inform the choice of data learning techniques for prediction of groundwater recharge in other places, while the obtained results will be instrumental for the development of future groundwater management strategies in the region.

Identification of important but unintuitive influential factors, such as Rainfall4-10, MeDS, MeWS, Mintem and Mintem5-9 can be used to inform monitoring requirements and improve the design of policy and management plans.

5.4. Advantages, Limitations and Further Research

These experiments for temporal prediction of groundwater recharge clearly demonstrate the effectiveness and robustness of the deep learning approach. The performance of the LSTM model appeared to be almost unaffected when various ratios of trained dataset were considered, or when influential predictors were removed.

Limitations of the study include: non-spatial data sources, short prediction time frames and quantification of uncertainty. Spatial data inputs have not been used in this paper to predict groundwater recharge. In future research, spatial data such as soil type, vegetation type, and slope may be used for improving the model predictions of groundwater recharge. Finally, our research does not yet quantify the uncertainty in prediction results. Although uncertainty analysis methods have been widely applied in classical prediction models [50,51,52], they have not generally been used in machine learning methods, especially deep learning models. Future work will introduce uncertainty analysis methods into these models for the uncertainty quantification of model predictions.

In future studies, nature-inspired intelligent algorithms [53,54] and other neural networks [55] will also be tested in order to further improve the fitness of predicted data. The performance of the models may also be enhanced by incorporating other environmental data (e.g., soil data) or supplemented interpolation data [56]. Furthermore, the iterative calculation [57] may also be used better optimize the model parameters. For this study, the whole observed groundwater recharge data set is divided into trained and predicted subsets based on various ratios. The trained dataset and predicted dataset were used for model training and validation respectively. In the future work, longer term predictions of groundwater recharge may be modelled, which can then be applied for to forecast recharge and managing regional groundwater more effectively.

6. Conclusions

The linear regression, multi-layer perception and long short-term memory models were used for the challenging problem of predicting long-term, time-continuous groundwater recharge. The LSTM model greatly outperformed the MLP and linear regression models when the ratio of the training dataset to the full dataset (composed of the training dataset and the prediction dataset) was set as 70%, 80% or 90%. In turn, the MLP model outperformed the linear regression model. The results clearly demonstrate the effectiveness of the three models for various ratios of trained data to predicted data. The relative importance of potential predictors associated with the observed variation in regional groundwater recharge was also assessed.