Groundwater Depth Prediction Using Data-Driven Models with the Assistance of Gamma Test

Jiyang Tian 1, Chuanzhe Li 1,*, Jia Liu 1,2, Fuliang Yu 1, Shuanghu Cheng 3, Nana Zhao 4 and Wan Zurina Wan Jaafar 5 1 State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China; tjyshd@126.com (J.T.); hettyliu@126.com (J.L.); yufl@iwhr.com (F.Y.) 2 State Key Laboratory of Hydrology-Water Resource and Hydraulirc Engineering, Hohai University, Nanjing 210098, China 3 Bureau of Water Resources Survey of Heibei, Shijiazhuang 050031, China; hbcsh@163.com 4 Institute of Wetland Research, Chinese Academy of Forestry, Beijing 100091, China; annazhao2009@163.com 5 Department of Civil Engineering, Faculty of Engineering, University of Malaya, Kuala Lumpur 50603, Malaysia; wzurina@um.edu.my * Correspondence: azhe051@163.com; Tel.: +86-10-6878-1656


Introduction
Groundwater is an important part of water resources for domestic, agricultural, industrial, and environmental uses due to its generally good quality, widespread availability, good water stability and the fact it is not easily polluted [1].Groundwater resources are especially important in the arid and semi-arid regions, which can cover the shortage of water caused by uneven distribution of the surface water in time and space.However, arable land increases, urban expansion and population growth have caused a dramatic increase in water consumption, which has led to a decrease of the groundwater level in many countries all over the world [2,3].
With the rapid economic development in China in the past several decades, the demand for water has increased at a fast speed.This has led to an excessive use of water resources, especially groundwater.At present, there are many groundwater overdraft areas in China, among which Hebei Province in North China plain is one of the most serious overdraft areas where the groundwater level has decreased dramatically [4,5].Much attention has been paid by both researchers and policy makers in this area.There have been studies focusing on both the reasons for and the negative effects of the continuous decline of the groundwater level in Hebei province and the North China plain [6,7].
Groundwater dynamic forecasting provides important information for groundwater management.Physical descriptive models and data-driven models are two classes of dynamic prediction models [8].Detailed hydrogeological data in space and time, which are usually obtained by a large number of field experiments, are required to construct the physical models in practical application.So, a data-driven model is a good choice when data are limited and the groundwater system is complex.Power Function Model (PFM), Back-Propagation Artificial Neural Network (BPANN) and Support Vector Machines (SVM) are popular data-driven models used for forecasting the groundwater level [9][10][11].Among these methods, PFM is a simple non-linear regression (NLR) model, while the other two models are more complex.Some studies indicate that the BPANN model and SVM model perform better, especially the SVM model [12], but the prediction abilities of these two models are not always good for all cases due to the complexities of hydrogeological conditions and the groundwater flow in different regions [13].
There are many causes for the decrease of the groundwater level, including natural factors, anthropic factors, biological factors and economic factors, which can all be treated as the inputs of the data-driven models.To construct a reliable data-driven model for groundwater prediction, some challenging questions need to be answered beforehand.For example, which inputs are relevant to the groundwater level and which are irrelevant?To what extent do the inputs determine the output from a smooth model?How many data points are required to make a prediction with the best possible accuracy?So far, these questions have not been addressed adequately in groundwater level prediction using data-driven models.However, it is possible that significant progress can be made to tackle these questions, because the Gamma Test (GT) has been presented as an efficient algorithm.By estimating the variance of the noise in the raw data, GT can directly estimate the best mean squared error for a given selection of input by a smooth model on unseen data, before the model construction.In this way, it can help find the best input data combination and identify the sufficient number of data points used for model training.A formal proof of the GT can be found in Evans and Jones [14].Then, a series of studies have indicated that GT can greatly reduce the model construction work and can potentially be used in the area of hydrological prediction and water management.Moghaddamnia and colleagues used GT to select the input data for Artificial Neural Networks (ANN) and the adaptive neuro-fuzzy inference system for evaporation estimation [15].Noori and colleagues assessed the SVM model performance for monthly stream flow prediction with the assistance of GT [16].Han and colleagues explored the model structure for index flood regionalization using GT [17].However, with respect to groundwater prediction, little work has been done using this efficient tool.
In this study, GT is used to rank the relative importance of the model inputs and find the best input combinations before the data-driven models are calibrated.Three data-driven models, i.e., PFM, BPANN and SVM, are used to predict the variations of the groundwater depth and the performances of the three established models are further compared.Twelve indices, including natural, anthropic, biological, economic and social factors which may influence the groundwater depth, are considered as the input of the data-driven models.The study is carried out in the plain of Shijiazhuang, the capital city of Hebei province in North China.During the past three decades, the groundwater level in this region has declined by 30 m, which has not only restricted the development of the economy but has also caused serious environmental problems [18].It is hoped that the study can provide a novel and complementary methodology for the prediction of groundwater level dynamics in the study area.

Study Area
Shijiazhuang plain (from lat 37 • 33 N to lat 38 • 42 N and from long 113 • 18 E to long 114 • 41 E) is located in central and eastern Hebei province and belongs to the Bohai Sea Economic Zone.The total area of Shijiazhuang plain is 8157 km 2 .The average annual precipitation is about 490 mm, and the distribution is quite uneven in space and time.Over 70% of the precipitation occurs in the summer months from June to August.The average air temperature ranges from −0.

Study Area
Shijiazhuang plain (from lat 37°33′ N to lat 38°42′ N and from long 113°18′ E to long 114°41′ E) is located in central and eastern Hebei province and belongs to the Bohai Sea Economic Zone.The total area of Shijiazhuang plain is 8157 km 2 .The average annual precipitation is about 490 mm, and the distribution is quite uneven in space and time.Over 70% of the precipitation occurs in the summer months from June to August.The average air temperature ranges from −0.8 °C in the winter seasons to 25.9 °C in the summer seasons.The plain comprises Shijiazhuang city and 13 counties, including Xingtang, Luquan, Yuanshi, Gaoyi, Xinle, Zhengding, Luancheng, Zhaoxian, Wuji, Gaocheng, Jinzhou, Xinji and Shenze.As a semi-arid region, water resources are critical to economic development and groundwater is the main source of water supply.The location of the study area and the monitoring wells and meteorological stations located in this area are shown in Figure 1

Input Indices and Data Source
Natural, anthropic, biological, economic and social factors which are closely related to the decline of the groundwater level are chosen as the key inputs of the data-driven models.These are 12 indices, including precipitation, evaporation, groundwater exploitation, industrial groundwater use, agricultural irrigation groundwater use, crop yield, population, urbanization rate, primary industry output, secondary industry output, tertiary industry output and gross domestic product (GDP).Precipitation and evaporation are the main natural factors which influence the groundwater recharge and loss.Groundwater exploitation is the groundwater fetching and use quantity, and it is the main anthropic factor which is the major reason for groundwater descending in the study area in recent decades [19], while industrial and agricultural irrigation groundwater uses are two major causes of groundwater exploitation.Crop yield is the main biological factor which reflects the crop

Input Indices and Data Source
Natural, anthropic, biological, economic and social factors which are closely related to the decline of the groundwater level are chosen as the key inputs of the data-driven models.These are 12 indices, including precipitation, evaporation, groundwater exploitation, industrial groundwater use, agricultural irrigation groundwater use, crop yield, population, urbanization rate, primary industry output, secondary industry output, tertiary industry output and gross domestic product (GDP).Precipitation and evaporation are the main natural factors which influence the groundwater recharge and loss.Groundwater exploitation is the groundwater fetching and use quantity, and it is the main anthropic factor which is the major reason for groundwater descending in the study area in recent decades [19], while industrial and agricultural irrigation groundwater uses are two major causes of groundwater exploitation.Crop yield is the main biological factor which reflects the crop water consumption, especially the groundwater consumption in the study area [20].GDP is a comprehensive index indicating the regional economic strength and the pressure on groundwater use, while primary industry output, secondary industry output and tertiary industry output can reflect different aspects of economic development.Primary industry includes agriculture, forestry, animal husbandry and fishery.Secondary industry includes mining, manufacturing and construction.Tertiary industry refers to all other economic activities not included in the primary or secondary industry.Population and urbanization rate are two important factors for social development.The growth of population and the increase of urbanization rate can also raise the use of water, especially groundwater in this region.Table 1 shows the serial number of the 12 indices and their mutual relations are illustrated by Figure 2.
Groundwater depth is the only index which is used to show the change of groundwater resources in this study.The higher the groundwater depth, the more scarce the groundwater resources.As located in Figure 1, 14 monitoring wells covering the whole study area are chosen to provide observational data on the annual groundwater depth from 1984 to 2013.The continuous decrease of groundwater depth can be easily observed over the 30-year period (Figure 3), and the seasonal variation of groundwater depth is obvious.The lowest groundwater level appears in summer, while the highest groundwater level can be found in winter.The groundwater depth is highly related to groundwater exploitation in different seasons for one year.In order to evaluate the groundwater of the whole study area, the average groundwater depth of 14 monitoring wells in each year is used by the models.Observed precipitation and evaporation at 13 meteorological stations are also obtained for the same period from the National Meteorological Information Center of China Meteorological Administration (available at http://www.nmic.gov.cn/).The average precipitation and evaporation from the 13 meteorological stations in each year are used as the model inputs.Thirty-year data of the other 10 input indices, including annual groundwater exploitation, industrial groundwater use, agricultural irrigation groundwater use, crop yield, population, urbanization rate, primary industry output, secondary industry output, tertiary industry output and GDP for the whole study area are also obtained for the period from 1984 to 2013.For each of the 10 input indices, annual value is the sum of the corresponding values of Shijiazhuang city and 13 counties.All the data of the 10 input indices come from the Shijiazhuang Statistical Yearbook [21].

Gamma Test
The GT estimates the mean square error (MSE) that can be achieved when modelling the unseen data using any continuous nonlinear model.A formal proof of GT is given by Evans and Jones [14].The basic idea is quite distinct from the earlier attempts with nonlinear analysis.Suppose we have a set of data observations in the form: where xi ∈ R M are input vectors confined to some closed bounded set C ∈ R M and, yi ∈ R are corresponding outputs.The system of GT can be expressed as the following form:

Gamma Test
The GT estimates the mean square error (MSE) that can be achieved when modelling the unseen data using any continuous nonlinear model.A formal proof of GT is given by Evans and Jones [14].The basic idea is quite distinct from the earlier attempts with nonlinear analysis.Suppose we have a set of data observations in the form: where xi ∈ R M are input vectors confined to some closed bounded set C ∈ R M and, yi ∈ R are corresponding outputs.The system of GT can be expressed as the following form:

Gamma Test
The GT estimates the mean square error (MSE) that can be achieved when modelling the unseen data using any continuous nonlinear model.A formal proof of GT is given by Evans and Jones [14].The basic idea is quite distinct from the earlier attempts with nonlinear analysis.Suppose we have a set of data observations in the form: where x i ∈ R M are input vectors confined to some closed bounded set C ∈ R M and, y i ∈ R are corresponding outputs.The system of GT can be expressed as the following form: where f is a smooth function and r is a random variable representing noise.Generally, the mean of the distribution of r is assumed as 0 and the variance of the noise Var(r) is bounded.The Gamma statistic Γ is the main parameter, which can estimate the model's output variance.
For each vector The distance statistic of input data can be calculated: where | . . .| denotes Euclidean distance, and the corresponding Gamma function of the output values: where y is the corresponding y-value for the kth nearest neighbor of x i in Equation ( 3).In order to compute Γ, the p points (δ M (k), γ M (k)) are calculated by univariate linear regression equation with least-squares: The value of Γ is the intercept of the Equation ( 5).Usually, the data, which are tested by GT, are normalized to a range of 0-1.Particularly, when δ equals 0, Γ value can be computed by γ: A more detailed description of calculation principles can be found in Evans and Jones [14].The V ratio is another term, which can return a scale invariant noise estimate: where σ 2 (y) is the variance of output y.According to the definition of V ratio , the value of V ratio close to 0 indicates a high degree of predictability of the given output y.In addition, the estimation of noise variance on the given output can be more credible if the standard error (SE) value is close to 0, and the complexity of the smooth function can be measured by the value of gradient.

Power Function Model (PFM)
The PFM is a kind of non-linear regression model, which is often used for prediction by establishing the relationship between the forecast factors and the forecast objects.The general equation is as follows: where y is the forecast object, p k (k = 0, . . .,N) are the parameters generally estimated by least squares and x k (k = 1, . . .,M) are the explanatory variables (forecast factors).

Back-Propagation Artificial Neural Network (BPANN)
The ANN is a nonlinear arithmetic system with densely interconnected processing elements or neurons.Three kinds of layers are contained in the mathematical structure, including input, hidden, and output layers.These layer has their nodes and activation functions.The back-propagation algorithm can effectively train the network and shorten the learning time, and it needs little information about the complex mechanism and process which should be explicitly described in mathematical form.If the neuron is the jth one in the present layer, while the inputs which it receives from the other n neurons are x 1 , x 2 , . . . . . .,x n , respectively, in the previous layer.The connection weights between the jth neuron and the other n neurons are w 1j , w 2j , . . . . . ., w nj , respectively.The mathematical expression is as follows: where y j is the output of the neuron, b j is the threshold of the neuron, and f is the transfer function.
In this study, one hidden layer is used to build the ANN model, which is trained by the back-propagation algorithm.A log-sigmoid and a linear function are contained in hidden layer and output layer, respectively.This kind of configuration is the most commonly used for ANNs, which can improve extrapolation ability.The input layer, hidden layer and output layer have five nodes, 30 nodes and one node, respectively (Figure 4).A momentum term is also added in the weight updating process to avoid the results being captured in local minimum.The number of hidden layer and nodes is discussed in Section 3.5.
between the jth neuron and the other n neurons are w1j, w2j, ……, wnj, respectively.The mathematical expression is as follows: where yj is the output of the neuron, bj is the threshold of the neuron, and f is the transfer function.
In this study, one hidden layer is used to build the ANN model, which is trained by the backpropagation algorithm.A log-sigmoid and a linear function are contained in hidden layer and output layer, respectively.This kind of configuration is the most commonly used for ANNs, which can improve extrapolation ability.The input layer, hidden layer and output layer have five nodes, 30 nodes and one node, respectively (Figure 4).A momentum term is also added in the weight updating process to avoid the results being captured in local minimum.The number of hidden layer and nodes is discussed in Section 3.5.

Support Vector Machines (SVMs)
The SVM model is based on Vapnik-Chervonenkis (VC) dimension and structural risk minimum principle [22].The SVM provides a new approach to solving the nonlinear and high dimension problem with a small sample set.The basic idea of SVM is to search for the nonlinear relationship between input and output vectors through nonlinear transformation of the input vector into the highdimensional feature space.Given a set of N samples of x is the input vector, k y is the corresponding output value), and the regression function of SVM can be expressed as: (10) where w is a weight vector, ϕ is a nonlinear transfer function that implements transformation the nonlinear to linear relationship of input to output vectors, and b is a bias.Vapnik introduced the

Support Vector Machines (SVMs)
The SVM model is based on Vapnik-Chervonenkis (VC) dimension and structural risk minimum principle [22].The SVM provides a new approach to solving the nonlinear and high dimension problem with a small sample set.The basic idea of SVM is to search for the nonlinear relationship between input and output vectors through nonlinear transformation of the input vector into the high-dimensional feature space.Given a set of N samples of {x k , y k } N k=1 (x k is the input vector, y k is the corresponding output value), and the regression function of SVM can be expressed as: where w is a weight vector, ϕ is a nonlinear transfer function that implements transformation the nonlinear to linear relationship of input to output vectors, and b is a bias.Vapnik introduced the convex quadratic optimization question to ensure that extreme solution is optimal [22], and a ε-insensitively loss function is added to Equation (10): where ξ and ξ * are slack variables that penalize training errors by the loss function over the error tolerance ε, and C denotes the degree of penalization for the sample outside the error tolerance ε.
The input vectors are called support vectors.The dual Lagrangian form is: with the constraints, where α and α* are Lagrange multipliers, and the optimal desired weight vector of the regression hyperplane can be expressed with the kernel function: where K(x, x i ) is the kernel function.Equation ( 10) can be expressed as: Linear kernel function (LKF), polynomial kernel function (PKF) and radial basis function (RBF) are commonly used as kernel functions.Among these three kernel functions, the number of hyper-parameters in PKF is much more than LKF and RBF, so LKF and RBF are chosen in this study [23,24].The parameters of RBF are discussed in Section 3.5.LKF, PKF and RBF are described by the Equations ( 16)- (18), respectively:

Implementation and Assessment of Three Models
In this study, we define y as the groundwater depth which is an index of the groundwater dynamic variations and the forecast object, and x k as the input factors which may influence the groundwater in the models of PFM, BPANN and SVM.PFM is constructed by the Statistical Product and Service Solutions (SPSS) software in version 19.0.BPANN and SVM are both implemented by MATLAB toolbox.
Twenty groups of data from 1984 to 2003 are used to build and train the models, while 10 groups of data from 2004 to 2013 are used for testing.All the data are normalized to a range of 0-1 to avoid disturbance of dimensions.The root mean squared error (RMSE) and Nash-Sutcliffe efficiency represented by E are used to assess the accuracy of the three built models based on the observed values and the simulated values: where y i and ŷi are the observed and simulated values of groundwater depth, respectively, y i is the mean observed value, i = 1,2, . . ., N, and N is the number of data groups.The perfect scores of RMSE and E are 0 and 1, respectively.

Correlation Test for 12 Inputs
Pearson correlation coefficients (ρ) are used to analyze the correlations among the 12 inputs.The value of 1 means a perfect positive correlation while the value of −1 means a perfect negative correlation.The correlation coefficients of the 12 inputs are shown in Table 2. Most correlation coefficients are lower than 0.5 and most inputs show weak correlations.However, the correlation coefficient between GDP and secondary industry output reach 0.799, which indicates that secondary industry output is the main means to improve the economy in Shijiazhuang plain.Tertiary industry output also has a high correlation with population, which means that tertiary industry is greatly influenced by population in the study area.Though the two pairs of inputs have relatively high correlation coefficients, overall, the correlations among the 12 inputs have little effect on the feature selection results.

Relative Importance of the Model Inputs
The GT can provide the least model error compared to any smooth model by only giving the input and output data.In the study, the Gamma value Γ is used as the main criterion for distinguishing the importance of model inputs and the other three factors produced by GT (i.e., Gradient, Standard error of Γ and V ratio ) are used as a reference.This is because there are still some technical problems of GT that cannot be explained in detail, such as the utilization of the three factors (Gradient, Standard error of Γ and V ratio ).Basically, the least |Γ| represents the best model input combination.The gradient can indicate the model complexity.The lower the gradient, the more simple the model should be fitted.The SE of Γ shows the reliability of the Gamma value.The higher it is, the more unreliable the Gamma value is.V ratio illustrates the predictability of outputs using the effective inputs [15].
A mathematical model (e.g., PFM, BPANN and SVM in this study) may be built in high quality, if the three factors (gradient, SE of Γ and V ratio ) show low values by setting the best input data [17].Table 3 shows the GT results for schemes with different combinations of model input indices listed in Table 1.Scheme 1 involves all 12 input indices and the other schemes have 11 indices with the rest masked.It is interesting to find that only Scheme 5 (without index 4, i.e., industrial groundwater use) and Scheme 13 (without index 12, i.e., urbanization rate) have larger |Γ| values than Scheme 1.
The other schemes all have smaller |Γ| values which indicate they will perform better than Scheme 1 in constructing a reliable model.It can be seen from Table 3 that the best combination of indices is in Scheme 3 which has the least |Γ| (0.00170), and the worst combination of indices is in Scheme 13 which has the highest |Γ| (0.01372).Then, we can rank the schemes from best to worst, which is 3 > 2 > 11 > 12 > 7 > 9 > 8 > 6 > 10 > 4 > 1 > 5 > 13, and the missing index in each of the schemes can form a ranking of the importance of the model input indices, which is 2 < 1 < 10 < 11 < 6 < 8 < 7 < 5 < 9 < 3 < 4 < 12.This indicates that natural factors (precipitation and evaporation) have less effect on the groundwater level, while anthropic factors, biological factors, economic factors and social factors have greater influence.Urbanization rate and industrial groundwater use are the two most important model input indices which imply that the process of modernization and industrial development in Shijiazhuang plain is the main cause of groundwater recession.

Model Input Selection Based on Gamma Test
In order to find the best combination of model inputs, search methods of backward and forward selections are coupled with the GT [25].The backward procedure starts with all indices and each index is gradually removed in turn, while the forward procedure starts with the minimum index number and carries on by adding one index at a time successively.The index which works the best with the indices from the last round (i.e., resulting in the least |Γ| value) is added in forward selection and the one with the absence of which leads to the best results is removed in backward selection.The procedure is iterated until one index is left for backward selection or until all indices are included for forward selection.Tables 4 and 5 show the model input selection results by using the search methods of backward and forward selections coupled with the GT, respectively.It can be found from the two tables that Scheme 19 has relatively fewer indices and the smallest values of |Γ| (0.00006).Theoretically, because of the least |Γ| and relatively fewer indices, Scheme 19 is the best model input combination.In this study, Scheme 19 is used to construct the three models.However, GT is still a method which requires more experiments to be completed.In order to make our conclusions more reliable, all the schemes in Tables 4 and 5 are calculated as part of the sensitivity analysis of the input choice in Section 3.6.

Influence of the Parameters in BPANN and SVM (RBF)
Possible data overfitting is a disadvantage of BPANN and SVM (RBF).In order to avoid this problem and make the models more accurate, we choose the hidden layers and nodes for BPANN and two parameters (C and γ) for SVM (RBF) carefully by means of the try-and-error method.The parameter C which is a positive trade-off parameter determines the degree of empirical error and the parameter γ which is the main parameter in kernel function of RBF.Numerous trials were conducted.
Taking Scheme 19 as an example, the optimum numbers of hidden layers and nodes per layer are 1 and 30 for BPANN, respectively, while the optimal values of C and γ in this study are chosen to be 1.0 and 0.3 for SVM (RBF), respectively.All these optimal parameters are chosen based on the values of RMSE and E. With the change in number of hidden layers and nodes, different values of RMSE and E can be seen in Figures 5 and 6.The influences of the C and γ on the testing results are shown in Figures 7 and 8.

Influence of the Parameters in BPANN and SVM (RBF)
Possible data overfitting is a disadvantage of BPANN and SVM (RBF).In order to avoid this problem and make the models more accurate, we choose the hidden layers and nodes for BPANN and two parameters (C and γ) for SVM (RBF) carefully by means of the try-and-error method.The parameter C which is a positive trade-off parameter determines the degree of empirical error and the parameter γ which is the main parameter in kernel function of RBF.Numerous trials were conducted.
Taking Scheme 19 as an example, the optimum numbers of hidden layers and nodes per layer are 1 and 30 for BPANN, respectively, while the optimal values of C and γ in this study are chosen to be 1.0 and 0.3 for SVM (RBF), respectively.All these optimal parameters are chosen based on the values of RMSE and E. With the change in number of hidden layers and nodes, different values of RMSE and E can be seen in Figures 5 and 6.The influences of the C and γ on the testing results are shown in Figures 7 and 8

Fitting Results and Model Errors
For Scheme 19, the equation of PFM can be expressed as follows calibrated on the 20 groups of data from 1984 to 2003, respectively: The other two models of BPANN and SVM are also trained based on the dataset from 1984 to 2003.The procedure is carried out several times for both models and the optimal results are adopted.The RMSE and E of the three models for fitting results are shown in Table 6.

Fitting Results and Model Errors
For Scheme 19, the equation of PFM can be expressed as follows calibrated on the 20 groups of data from 1984 to 2003, respectively: The other two models of BPANN and SVM are also trained based on the dataset from 1984 to 2003.The procedure is carried out several times for both models and the optimal results are adopted.The RMSE and E of the three models for fitting results are shown in Table 6.

Fitting Results and Model Errors
For Scheme 19, the equation of PFM can be expressed as follows calibrated on the 20 groups of data from 1984 to 2003, respectively: The other two models of BPANN and SVM are also trained based on the dataset from 1984 to 2003.The procedure is carried out several times for both models and the optimal results are adopted.The RMSE and E of the three models for fitting results are shown in Table 6.where x 3 is groundwater exploitation, x 4 is industrial groundwater use, x 6 is crop yield, x 9 is secondary industry output and x 12 is urbanization rate, ŷ is the fitting results for groundwater depth.
The other two models of BPANN and SVM are also trained based on the dataset from 1984 to 2003.The procedure is carried out several times for both models and the optimal results are adopted.The RMSE and E of the three models for fitting results are shown in Table 6.The RMSEs of PFM and SVM (LKF) are both higher than 1.0, and the value of SVM (LKF) is close to 1.5.The RMSE is lower than 0.5 for the SVM (RBF) model, while the value of BPANN model is higher than 0.5 and lower than PFM and SVM (LKF).So, the fitting values of SVM (RBF) are closer to observed values.The highest values of E are near 1.0 for SVM (RBF), and the lowest values of E are lower than 0.85 for SVM (LKF).It indicates that SVM (RBF) has the highest reliability in the four models.
Based on the values of RMSE and E, it can be easily observed that the SVM (RBF) model has the best generalization ability in the sample learning process, while the SVM (LKF) model performs the worst.Good fitting results do not mean good prediction ability for the data-driven models.The performance of the models needs to be tested with further datasets.The observed and fitted values (Scheme 19) of groundwater depth are shown in Figure 9.The RMSEs of PFM and SVM (LKF) are both higher than 1.0, and the value of SVM (LKF) is close to 1.5.The RMSE is lower than 0.5 for the SVM (RBF) model, while the value of BPANN model is higher than 0.5 and lower than PFM and SVM (LKF).So, the fitting values of SVM (RBF) are closer to observed values.The highest values of E are near 1.0 for SVM (RBF), and the lowest values of E are lower than 0.85 for SVM (LKF).It indicates that SVM (RBF) has the highest reliability in the four models.Based on the values of RMSE and E, it can be easily observed that the SVM (RBF) model has the best generalization ability in the sample learning process, while the SVM (LKF) model performs the worst.Good fitting results do not mean good prediction ability for the data-driven models.The performance of the models needs to be tested with further datasets.The observed and fitted values (Scheme 19) of groundwater depth are shown in Figure 9.

Testing Results and Model Errors
The three types of data-driven models are tested based on the 10 groups of data from 2004 to 2013 for three schemes.The RMSE and E of the three models for testing results are shown in Table 7.For Scheme 19, only the RMSE of SVM (RBF) is lower than 2.0, and the RMSE of SVM (LKF) is higher than 3.0.The BPANN model has lower RMSE than the PFM model.The value of E is higher than 0 for SVM (RBF), and the lowest values of E are lower than −2.0 for PFM and SVM (LKF).It can be noticed that BPANN and SVM (RBF) models perform better than PFM and SVM (LKF) as a whole with a relatively lower RMSE and better E, but stability is poor for the BPANN model.The SVM (RBF) model performs best in the testing process, and it has generalization ability in the testing process.The observed and tested values (Scheme 19) of groundwater depth are shown in Figure 10.

Testing Results and Model Errors
The three types of data-driven models are tested based on the 10 groups of data from 2004 to 2013 for three schemes.The RMSE and E of the three models for testing results are shown in Table 7.For Scheme 19, only the RMSE of SVM (RBF) is lower than 2.0, and the RMSE of SVM (LKF) is higher than 3.0.The BPANN model has lower RMSE than the PFM model.The value of E is higher than 0 for SVM (RBF), and the lowest values of E are lower than −2.0 for PFM and SVM (LKF).It can be noticed that BPANN and SVM (RBF) models perform better than PFM and SVM (LKF) as a whole with a relatively lower RMSE and better E, but stability is poor for the BPANN model.The SVM (RBF) model performs best in the testing process, and it has generalization ability in the testing process.The observed and tested values (Scheme 19) of groundwater depth are shown in Figure 10.

Sensitivity Analysis of the Input Choice
In order to make the input choice conclusion more reliable, the same methods are used to calculate the values of RMSE and E for other schemes in Tables 4 and 5, and the results are compared with Scheme 19.The optimal models for BPANN and SVM (RBF) are also chosen by testing the influence of the parameters in the models.The fitting and testing results and errors of PFM, BPANN and SVM models are shown in Tables 8 and 9. Figures 11 and 12 show the changes of RMSE and E with |Γ|.It indicates that Scheme 19 is the best model input combination with lower values of RMSE and E for the three models, and it also has relatively fewer indices.Overall, with the increase of the |Γ|, the performances of the three models become worse, so GT can be used to confirm the best combination of model inputs.It should be mentioned that Schemes 21, 22, 23, 24 and 26 with much fewer indices have much higher RMSE and lower E. So, the schemes with much fewer indices but smaller values of |Γ| should be tested carefully.Additionally, it also can be seen that the SVM (RBF) performs best in most schemes.

Sensitivity Analysis of the Input Choice
In order to make the input choice conclusion more reliable, the same methods are used to calculate the values of RMSE and E for other schemes in Tables 4 and 5, and the results are compared with Scheme 19.The optimal models for BPANN and SVM (RBF) are also chosen by testing the influence of the parameters in the models.The fitting and testing results and errors of PFM, BPANN and SVM models are shown in Tables 8 and 9. Figures 11 and 12 show the changes of RMSE and E with |Γ|.It indicates that Scheme 19 is the best model input combination with lower values of RMSE and E for the three models, and it also has relatively fewer indices.Overall, with the increase of the |Γ|, the performances of the three models become worse, so GT can be used to confirm the best combination of model inputs.It should be mentioned that Schemes 21, 22, 23, 24 and 26 with much fewer indices have much higher RMSE and lower E. So, the schemes with much fewer indices but smaller values of |Γ| should be tested carefully.Additionally, it also can be seen that the SVM (RBF) performs best in most schemes.

Discussion and Conclusions
Considering limited data and complicated geological conditions in the plain of Shijiazhuang, physical models are not a good choice, and three types of data-driven models are used in this study to predict the groundwater depth in Shijiazhuang plain.Five kinds of factors (natural, anthropic, biological, economic and social) including 12 indices which may influence the groundwater depth are considered, i.e., precipitation, evaporation, groundwater exploitation, industrial groundwater use, agricultural irrigation groundwater use, crop yield, population, urbanization rate, primary industry output, secondary industry output, tertiary industry output and GDP.An efficient tool named GT is adopted in this study to identify the relative importance of the indices and further determine the best

Discussion and Conclusions
Considering limited data and complicated geological conditions in the plain of Shijiazhuang, physical models are not a good choice, and three types of data-driven models are used in this study to predict the groundwater depth in Shijiazhuang plain.Five kinds of factors (natural, anthropic, biological, economic and social) including 12 indices which may influence the groundwater depth are considered, i.e., precipitation, evaporation, groundwater exploitation, industrial groundwater use, agricultural irrigation groundwater use, crop yield, population, urbanization rate, primary industry output, secondary industry output, tertiary industry output and GDP.An efficient tool named GT is adopted in this study to identify the relative importance of the indices and further determine the best

Discussion and Conclusions
Considering limited data and complicated geological conditions in the plain of Shijiazhuang, physical models are not a good choice, and three types of data-driven models are used in this study to predict the groundwater depth in Shijiazhuang plain.Five kinds of factors (natural, anthropic, biological, economic and social) including 12 indices which may influence the groundwater depth are considered, i.e., precipitation, evaporation, groundwater exploitation, industrial groundwater use, agricultural irrigation groundwater use, crop yield, population, urbanization rate, primary industry output, secondary industry output, tertiary industry output and GDP.An efficient tool named GT is adopted in this study to identify the relative importance of the indices and further determine the best model input combination.The results of GT indicate that natural factors (precipitation and evaporation) have less effect on the groundwater level, while anthropic factors, biological factors, economic factors and social factors are the main reasons for the groundwater level decline, especially the anthropic factors which reflect the groundwater exploitation.Evaporation is found to be the least relevant index among the 12 indices, while urbanization rate and industrial groundwater use are shown to have the most significant effect on groundwater variations.This indicates that the relationship between the evaporation and the groundwater becomes weaker with the drop in groundwater level.It also reveals that modernization and industrial development in Shijiazhuang plain are the main reasons for groundwater recession.
The RMSE and E are used to assess the accuracy of the data-driven models.The results show that BPANN and SVM (RBF) models perform better than PFM and SVM (LKF) as a whole in all three schemes, while SVM (RBF) has the best generalization ability during the fitting and testing process.The SVM (RBF) model performs the best in both model fitting and model testing.Additionally, all the schemes in Tables 4 and 5 are used to analyze the sensitivity of the input choice.The results show that GT can be used to confirm the best combination of model inputs, but schemes with much fewer indices and smaller values of |Γ| should be tested carefully.
The groundwater level has decreased dramatically and is continuing to do so in Shijiazhuang plain.Countermeasures need to be taken to reduce groundwater usage and guarantee the sustainable utilization of groundwater, by improving the water use efficiency, promoting comprehensive water-saving measures, adjusting the industrial structure and finding alternative water sources, such as water recycling and the South-to-North Water Diversion Project.The results of this study may also provide a methodology and a reference for the prediction of groundwater resources in data-limited areas, especially for the study region and North China in general.However, monthly prediction of the groundwater depth has not been considered in this study due to limited data observations.The advantage of GT can be better exploited when dealing with larger datasets using data-driven models in groundwater prediction.
8 • C in the winter seasons to 25.9 • C in the summer seasons.The plain comprises Shijiazhuang city and 13 counties, including Xingtang, Luquan, Yuanshi, Gaoyi, Xinle, Zhengding, Luancheng, Zhaoxian, Wuji, Gaocheng, Jinzhou, Xinji and Shenze.As a semi-arid region, water resources are critical to economic development and groundwater is the main source of water supply.The location of the study area and the monitoring wells and meteorological stations located in this area are shown in Figure 1.Two aquifers (the Holocene (Q4) and Upper Pleistocene (Q3)) in Shijiazhuang plain are composed of quaternary unconsolidated sediments, including sandy loam, sand and gravel-cobble.The west of Xingtang, Luquan and Yuanshi belongs to Q4, and other parts of the study area belong to Q3.The burial depth of the floor is 7-20 m for Q4 and 40-100 m for Q3.The specific field is 20-30 m 3 /(h•m) for Q4 and about 35 m 3 /(h•m) for Q3.Both of the aquifers contain a groundwater mineralization level lower than 0.5 g/L.All the 14 wells are used to monitor the shallow groundwater.
. Two aquifers (the Holocene (Q4) and Upper Pleistocene (Q3)) in Shijiazhuang plain are composed of quaternary unconsolidated sediments, including sandy loam, sand and gravel-cobble.The west of Xingtang, Luquan and Yuanshi belongs to Q4, and other parts of the study area belong to Q3.The burial depth of the floor is 7-20 m for Q4 and 40-100 m for Q3.The specific field is 20-30 m 3 /(h•m) for Q4 and about 35 m 3 /(h•m) for Q3.Both of the aquifers contain a groundwater mineralization level lower than 0.5 g/L.All the 14 wells are used to monitor the shallow groundwater.

Figure 1 .
Figure 1.The location of the study area and the distributions of monitoring wells and meteorological stations.

Figure 1 .
Figure 1.The location of the study area and the distributions of monitoring wells and meteorological stations.

Figure 2 .
Figure 2. Mutual relations of the 12 input indices.The blue box shows the five key factors, which are closely related to groundwater level changes.The green, red, cyan, grey and yellow boxes respectively represent natural, anthropic, economic, biological and social indices.The connecting lines with twoway arrows mean interaction effects between the indices and the groundwater level.The connecting lines with single arrows mean one-way supplement or consumption of the groundwater, and the lines without arrow mean inclusiveness.

Figure 3 .
Figure 3. Precipitation and groundwater depth from 1984 to 2013 in Shijiazhuang plain.

Figure 2 .
Figure 2. Mutual relations of the 12 input indices.The blue box shows the five key factors, which are closely related to groundwater level changes.The green, red, cyan, grey and yellow boxes respectively represent natural, anthropic, economic, biological and social indices.The connecting lines with two-way arrows mean interaction effects between the indices and the groundwater level.The connecting lines with single arrows mean one-way supplement or consumption of the groundwater, and the lines without arrow mean inclusiveness.

Figure 2 .
Figure 2. Mutual relations of the 12 input indices.The blue box shows the five key factors, which are closely related to groundwater level changes.The green, red, cyan, grey and yellow boxes respectively represent natural, anthropic, economic, biological and social indices.The connecting lines with twoway arrows mean interaction effects between the indices and the groundwater level.The connecting lines with single arrows mean one-way supplement or consumption of the groundwater, and the lines without arrow mean inclusiveness.

Figure 3 .
Figure 3. Precipitation and groundwater depth from 1984 to 2013 in Shijiazhuang plain.

Figure 3 .
Figure 3. Precipitation and groundwater depth from 1984 to 2013 in Shijiazhuang plain.

Figure 4 .
Figure 4.The architecture of the Back-Propagation Artificial Neural Network (BPANN) model with 12 input layers, 30 hidden layers and one output layer.

Figure 4 .
Figure 4.The architecture of the Back-Propagation Artificial Neural Network (BPANN) model with 12 input layers, 30 hidden layers and one output layer. .

Figure 5 .Figure 5 .
Figure 5. Influence of number of hidden layers on the BPANN model with 30 nodes.(a) RMSE; (b) E.

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Influence of number of nodes in hidden layer on the BPANN model with one hidden layer.(a) RMSE; (b) E.
x3 is groundwater exploitation, x4 is industrial groundwater use, x6 is crop yield, x9 is secondary industry output and x12 is urbanization rate, ŷ is the fitting results for groundwater depth.

Figure 6 .Figure 6 .Figure 7 .Figure 8 .
Figure 6.Influence of number of nodes in hidden layer on the BPANN model with one hidden layer.(a) RMSE; (b) E. (a) (b) Figure 6.Influence of number of nodes in hidden layer on the BPANN model with one hidden layer.(a) RMSE; (b) E.
x3 is groundwater exploitation, x4 is industrial groundwater use, x6 is crop yield, x9 is secondary industry output and x12 is urbanization rate, ŷ is the fitting results for groundwater depth.

Figure 6 .Figure 7 .Figure 8 .
Figure 6.Influence of number of nodes in hidden layer on the BPANN model with one hidden layer.(a) RMSE; (b) E.
x3 is groundwater exploitation, x4 is industrial groundwater use, x6 is crop yield, x9 is secondary industry output and x12 is urbanization rate, ŷ is the fitting results for groundwater depth.

3. 5 .
Fitting Results and Model ErrorsFor Scheme 19, the equation of PFM can be expressed as follows calibrated on the 20 groups of data from 1984 to 2003, respectively:

Figure 9 .
Figure 9.The observed and fitted values (Scheme 19) of groundwater depth from 1984 to 2003.

Figure 11 .Figure 12 .
Figure 11.The different values of RMSE and E with the change of |Γ| in fitting results.(a) RMSE; (b) E.

Figure 11 .
Figure 11.The different values of RMSE and E with the change of |Γ| in fitting results.(a) RMSE; (b) E.

Figure 11 .Figure 12 .
Figure 11.The different values of RMSE and E with the change of |Γ| in fitting results.(a) RMSE; (b) E.

Figure 12 .
Figure 12.The different values of RMSE and E with the change of |Γ| in testing results.(a) RMSE; (b) E.

Table 1 .
Serial number of the 12 input indices.

Table 3 .
GT results for schemes with different combinations of model input indices.

Table 4 .
Backward selection results with the assistance of GT.

Table 5 .
Forward selection results with the assistance of GT.

Table 5 .
Forward selection results with the assistance of GT.

Table 6 .
Fitting results and errors of PFM, BPANN and SVM models for Scheme 19.

Table 6 .
Fitting results and errors of PFM, BPANN and SVM models for Scheme 19.

Table 7 .
Testing results and errors of PFM, BPANN and SVM models for Scheme 19.
Figure 9.The observed and fitted values (Scheme 19) of groundwater depth from 1984 to 2003.

Table 7 .
Testing results and errors of PFM, BPANN and SVM models for Scheme 19.The observed and tested values (Scheme 19) of groundwater depth from 2004 to 2013.

Table 8 .
Fitting results and errors of PFM, BPANN and SVM models.

Table 9 .
Testing results and errors of PFM, BPANN and SVM models.The observed and tested values (Scheme 19) of groundwater depth from 2004 to 2013.

Table 8 .
Fitting results and errors of PFM, BPANN and SVM models.

Table 9 .
Testing results and errors of PFM, BPANN and SVM models.