A Comparative Study of Groundwater Level Forecasting Using Data-Driven Models Based on Ensemble Empirical Mode Decomposition

Yicheng Gong 1,2,*, Zhongjing Wang 1,2,3,*, Guoyin Xu 1 and Zixiong Zhang 1 1 Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China; xgy13@mails.tsinghua.edu.cn (G.X.); zx-zhang15@mails.tsinghua.edu.cn (Z.Z.) 2 State Key Laboratory of Hydro-Science and Engineering, Tsinghua University, Beijing 100084, China 3 State Key Lab of Plateau Ecology and Agriculture, Qinghai University, Xining 810016, China * Correspondence: gongyc@tsinghua.edu.cn (Y.G.); zj.wang@tsinghua.edu.cn (Z.W.); Tel.: +86-010-6278-2021 (Z.W.)


Introduction
Groundwater is an increasingly important water resource for irrigation, and domestic and industrial activities in many countries. More reliable and accurate estimation of groundwater levels can help prevent groundwater overexploitation and improve water-use efficiency for water resource management. However, in some regions, the groundwater has been pumped out much faster than it can be replenished, which eventually reduces the groundwater level. In addition, groundwater level time series are highly non-linear and non-stationary in nature, and prediction depends on many complex environmental factors, such as groundwater aquifers, precipitation, etc. Groundwater aquifers are intrinsically heterogeneous systems that are affected by complex hydrogeological conditions with groundwater-surface water interactions at various temporal-spatial scales [1,2]. Therefore, it is essential to develop more effective models for groundwater level prediction. Many groundwater modeling approaches and data-driven models have been implemented to forecast groundwater levels [2][3][4]. A groundwater numerical model is firstly developed from a conceptual model, which often includes only the main and fundamental principles but neglects the complexities of groundwater systems [1]. It is hard to establish groundwater flow equations and set the values of hydrogeological parameters for conceptual models. Also, it is difficult to obtain long time series data for groundwater numerical modeling; there are still challenges and uncertainties in modelling processes [1,5].
Recently, data-driven models such as artificial neural networks (ANN), support vector machines (SVM), adaptive neuro-fuzzy inference systems (ANFIS) and genetic programming (GP), and time series methods such as autoregressive moving average (ARMA) and autoregressive integrated moving average (ARIMA) have been proved efficient in forecasting hydrologic time series (e.g., groundwater level, water demand and inflow) [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Yoon et al. [5] developed ANN and SVM models to forecast groundwater level fluctuations in a coastal aquifer. The performance of the SVM was similar to or even better than that of the ANN model, according to the results from the validation. Shirmohammadi et al. [10] applied time series, system identification and ANFIS models to predict groundwater levels. The obtained results showed that ANFIS is superior to time series and system identification. Moosavi [11] proposed a model to forecast groundwater levels for different prediction periods. The achieved results showed that predictive ability of the wavelet-ANIFIS model was more accurate than those of the ANN, ANFIS and wavelet-ANN models. Shiri and Kişi [12] investigated the ability of adaptive neuro-fuzzy inference systems (ANFIS) and genetic programming (GP) to predict water table depth fluctuations. The results showed that the ANFIS and GP models can be applied successfully in groundwater depth prediction. Fallah-Mehdipour et al. [6] investigated the capability of ANFIS and GP to forecast and simulate groundwater levels in the Karaj plain of Iran. These results indicated that GP is an effective method for predicting groundwater levels. Valipour et al. [13] compared autoregressive moving average (ARMA) with autoregressive integrated moving average (ARIMA) models in forecasting the inflow of the Dez dam reservoir. The ARIMA model was proved have better predictive ability compared to the ARMA model. Xu et al. [14] built the instance-based weighting and support vector regression data-driven models to reduce the predictive error of groundwater models. The results of two real-world case studies showed that data-driven models can be applied effectively to reduce the root mean square error of the groundwater models. Asefa et al. [15] used support vector machines (SVMs) to monitor network design. The obtained results showed that SVMs can select the best configurations of well networks by reproducing the behavior of Monte Carlo flow.
Ensemble empirical mode decomposition (EEMD) is a scale-adaptive method proposed by Wu and Huang [19], which improved on empirical mode decomposition (EMD) to avoid the drawback of mode mixing [20]. EEMD has a high ability to decompose the original signal into intrinsic mode function (IMF) components and residual components for nonstationary and nonlinear signal sequences. In the hydrological and environmental research field, EEMD has been successfully applied to predict nonlinear problems such as runoff [21][22][23], wind speed [24], wave height [25], particulate matter 2.5 (PM2.5) [26], streamflow [27,28], vegetation dynamics [29], etc. Wang et al. [21] proposed an EEMD-ARIMA model for the forecasting of annual runoff time series. Quantitative standard statistical performance values proved that EEMD-ARIMA is a superior model to ARIMA for the forecasting of annual runoff. Liu et al. [24] applied the fast ensemble empirical mode decomposition-multilayer perceptron (FEEMD-MLP), Fast Ensemble Empirical Mode Decomposition-Adaptive neuro fuzzy inference system (FEEMD-ANFIS), wavelet packet-multilayer perceptron (WP-MLP) and WP-ANFIS models to forecast wind speed. The results suggested that all of the proposed hybrid algorithms are suitable for wind speed prediction. Duan et al. [25] developed the improved empirical model decomposition-support vector regression (EMD-SVR) model for the short-term prediction of wave height and found that the EMD-SVR model performs better than the wavelet-decomposition-based SVR (WD-SVR) model. Ausati and Amanollahi [26] used ensemble empirical mode decomposition-general regression neural network (EEMD-GRNN), principal component regression (PCR), adaptive neuro-fuzzy inference system (ANFIS) and multiple liner regression (MLR) models to predict concentrations of particulate matter 2.5 (PM2.5) in the city of Sanandaj. It was concluded that the EEMD-GRNN hybrid model had higher accuracy than the linear model in forecasting the particulate matter 2.5 (PM2.5) concentration. Karthikeyan and Nagesh [28] tested the predictability of monthly streamflow time series using wavelet-and EMD-based ARMA models. The result showed that the wavelet-based ARMA method outperformed the EMD-based ARMA method. Hawinkel et al. [29] used EEMD to decompose the time series of the normalized difference vegetation index (NDVI) and extracted climate-driven interannual vegetation dynamics. The result indicated that the EEMD was a feasible method for the assessment of climate-driven interannual vegetation dynamics.
ANN, SVN and ANFIS data-driven models and EEMD were applied in formal hydrology studies. All of these models were used alone or in combination to predict nonlinear and nonstationary time series. In this study, the ANN, SVN and ANFIS nonlinear time-series model coupled with EEMD were applied to groundwater level prediction, for validity and accuracy testing. The main advantage of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS hybrid models is that exogenous factors affecting the groundwater level do not need to be considered in the prediction. In order to compare the forecasting performance of the three hybrid models, single ANN, SVM and ANFIS models were applied to predict the groundwater level by considering exogenous factors, such as precipitation, temperature and surface water level, etc. The prediction results obtained from models with and without coupling were compared based on standard statistical analysis. This paper is organized as follows: Section 2 briefly describes the basic theory of EEMD, ANN, SVM and ANFIS, and presents the framework of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS hybrid models; Section 3 introduced the study area, available data and performance criteria. Section 4 explained the groundwater level decomposition by the EEMD for modeling and application; Section 5 analyzed the predicting results of the proposed EEMD-ANN, EEMD-SVM and EEMD-ANFIS hybrid models; Section 6 concluded this study.

Ensemble Empirical Mode Decomposition
Ensemble empirical mode decomposition (EEMD) is a new empirical mode decomposition (EMD)-based algorithm, which is a fully data-adaptive method to analyze nonlinear and nonstationary signals [19,20]. EMD can decompose the original signal into the number of intrinsic mode functions (IMFs) when the signal meets two conditions: (1) the functions are symmetric and the mean value of the upper and lower envelopes should be zero; (2) the number of extrema and the number of zero crossings must be equal or differ at most by one. However, EMD still has some drawbacks, such as mode mixing. Ensemble empirical mode decomposition (EEMD) was proposed to combine the white-noise-assisted system based on EMD to solve the problem of mode mixing [19,30]. If x(t) is the signal or time series to be decomposed, the EMD algorithm can be briefly summarized as follows: Step 1. Identify all local maxima and minima points of the time series x(t); Step 2. Interpolate between all local maxima and minima points of x(t) to form the upper envelope e max (t) and lower envelope e min (t); Step 3. Compute the mean envelope m(t) between the upper envelope e max (t) and the lower envelope e min (t); m(t) = (e max (t) + e min (t))/2 Step 4. Calculate the IMF candidate; Step 5. Determine whether or not h(t + 1) satisfies the two conditions of IMF. Is h(t + 1) an IMF?
If yes, save h(t + 1), calculate the residue r(t h(i), t = t + 1, and treat r(t) as input data in step 2. If no, treat h(t + 1) and as input data in step 2 until it satisfies the two conditions of IMF.
Step 6. Continue until the final residue meets some predefined stopping criteria.
At the end of the shifting process, the final residual term r n (t) has less than two local extrema, and the original signal or time series can be decomposed into a set of IMFs and a residue as follows: where n is the number of IMFs, r n (t) is the final residue and h i (t) as IMFs are nearly orthogonal to each other. Ensemble empirical mode decomposition (EEMD) was proposed by adding a finite amount of Gaussian white noise to series with EMD mixture, and the mode mixing is mostly eliminated. Based on the EMD method, the EEMD can be briefly summarized as follows [19,31]: Step 1. Initialize the ensemble number and the amplitude of the added white noise.
Step 2. Add random white noise to produce the noise-added data.
Step 3. Identify the local maxima and minima and obtain the upper and lower envelopes.
Step 4. Compute the mean of the upper and lower envelopes.
Step 5. Decompose the data with added random white noise into IMFs.
Step 6. Repeat step 3 to step 5 until the stopping criteria. After the shift processing, the IMFs and the residue are obtained.

Artificial Neural Network
An artificial neural network (ANN) is a black box tool that has certain performance characteristics that resemble the biological neural networks in the human brain [32]. Feed-forward neural networks (FFNN) are the most commonly used artificial neural networks, and have been broadly employed for modeling and forecasting in hydrogeology [33,34]. The architecture of the FFNN model consists of an input layer, hidden layer, output layer and artificial neurons (Figure 1). In a feed forward network, the neurons in each layer are connected to those in the next layer. Therefore, the output of a node in a layer only depends on the input it received from previous layer, determined weight, and type of transform function. The present research used the three-layer FFNN model, which was trained with a Bayesian regularization (BR) algorithm. The tan-sigmoid transfer function was selected by previous studies [35]. The mathematical expression can be expressed as follows: where x i is the input vector, y j is the output, b i is the bias, w ji is a weight connecting x i and y j , N is the number of nodes, and f is the activation function.

Support Vector Machine
A support vector machine (SVM) is a novel machine learning technique based on statistical learning theory [36,37]. In a regression SVM model, a regression hyperplane with a ε-insensitivity loss function is a convex dual optimization problem. The solution can be obtained from the optimization algorithm. The deterministic function of the SVM can be expressed as follows: where w i is a weight vector, and b is a bias. x is mapped to a high-dimensional feature space by nonlinear transfer function φ. In this study, the sequential minimal optimization (SMO) algorithm was used to solve the dual optimization problem of SVM [38,39]. The SMO algorithm can be directly acquired from an analytical solution without invoking a quadratic optimizer. The programming codes of the library for support vector machines (LIBSVM) was applied to predict a non-linear time series in the validation and calibration [40]. Figure 2 shows the schematic diagram of the SVM model.
where i w is a weight vector, and b is a bias. x is mapped to a high-dimensional feature space by nonlinear transfer function φ . In this study, the sequential minimal optimization (SMO) algorithm was used to solve the dual optimization problem of SVM [38,39]. The SMO algorithm can be directly acquired from an analytical solution without invoking a quadratic optimizer. The programming codes of the library for support vector machines (LIBSVM) was applied to predict a non-linear time series in the validation and calibration [40]. Figure 2 shows the schematic diagram of the SVM model.

Adaptive Neuro Fuzzy Inference System
An adaptive neuro fuzzy inference system (ANFIS) as a hybrid algorithm is an adaptive neural network based on a fuzzy inference system [41]. ANFIS is capable of approximating any real continuous function in a compact set to any degree of accuracy. This study used the first order Sugeno-fuzzy model of ANFIS [41,42]. To simplify, it was assumed that the fuzzy inference system has two inputs x and y and one output f . For the first-order Sugeno inference system, the typical rules can be expressed as: Rule 2: If x is A2 and y is B2; then f2 = p2x + q2y + r2 where x and y are the crisp inputs to the node i , i A and i B are the linguistic labels (low, medium, high, etc.,) characterized by the convenient membership functions. i p , i q , and i r ( ) 1, 2 i = are the parameters in the then part (consequent part) of the first-order Sugeno fuzzy model. The Sugeno-fuzzy structure of ANFIS comprised a five-layer network ( Figure 3). More comprehensive information about ANFIS can be found in the literature [11,[43][44][45].

Input Layer
Hidden Layer Output Layer where i w is a weight vector, and b is a bias. x is mapped to a high-dimensional feature space by nonlinear transfer function φ . In this study, the sequential minimal optimization (SMO) algorithm was used to solve the dual optimization problem of SVM [38,39]. The SMO algorithm can be directly acquired from an analytical solution without invoking a quadratic optimizer. The programming codes of the library for support vector machines (LIBSVM) was applied to predict a non-linear time series in the validation and calibration [40]. Figure 2 shows the schematic diagram of the SVM model.

Adaptive Neuro Fuzzy Inference System
An adaptive neuro fuzzy inference system (ANFIS) as a hybrid algorithm is an adaptive neural network based on a fuzzy inference system [41]. ANFIS is capable of approximating any real continuous function in a compact set to any degree of accuracy. This study used the first order Sugeno-fuzzy model of ANFIS [41,42]. To simplify, it was assumed that the fuzzy inference system has two inputs x and y and one output f . For the first-order Sugeno inference system, the typical rules can be expressed as: Rule 2: If x is A2 and y is B2; then f2 = p2x + q2y + r2 where x and y are the crisp inputs to the node i , i A and i B are the linguistic labels (low, medium, high, etc.,) characterized by the convenient membership functions. i p , i q , and i r ( ) 1, 2 i = are the parameters in the then part (consequent part) of the first-order Sugeno fuzzy model. The Sugeno-fuzzy structure of ANFIS comprised a five-layer network ( Figure 3). More comprehensive information about ANFIS can be found in the literature [11,[43][44][45].

Input Layer
Hidden Layer Output Layer Figure 2. Schematic diagram of a SVM.

Adaptive Neuro Fuzzy Inference System
An adaptive neuro fuzzy inference system (ANFIS) as a hybrid algorithm is an adaptive neural network based on a fuzzy inference system [41]. ANFIS is capable of approximating any real continuous function in a compact set to any degree of accuracy. This study used the first order Sugeno-fuzzy model of ANFIS [41,42]. To simplify, it was assumed that the fuzzy inference system has two inputs x and y and one output f . For the first-order Sugeno inference system, the typical rules can be expressed as: Rule 1: If x is A 1 and y is B 1 ; then f 1 = p 1 x + q 1 y + r 1 (6) Rule 2: If x is A 2 and y is B 2 ; then f 2 = p 2 x + q 2 y + r 2 where x and y are the crisp inputs to the node i, A i and B i are the linguistic labels (low, medium, high, etc.,) characterized by the convenient membership functions. p i , q i , and r i (i = 1, 2) are the parameters in the then part (consequent part) of the first-order Sugeno fuzzy model. The Sugeno-fuzzy structure of ANFIS comprised a five-layer network ( Figure 3). More comprehensive information about ANFIS can be found in the literature [11,[43][44][45].
Water 2018, 10, x FOR PEER REVIEW 6 of 20

The Hybrid EEMD-ANN, EEMD-SVM and EEMD-ANFIS Forecasting Models
The goal of this study was to improve the forecasting accuracy of groundwater levels by coupling the EEMD and data-driven models. A schematic diagram of the hybrid EEMD-ANN, EEMD-SVM and EEMD-ANFIS models is illustrated in Figure 4. As shown in Figure 4, the process of running the three hybrid models consists of three main steps:  (2) Secondly, the three data driven models ANN, SVM and ANFIS are developed to make the corresponding prediction using extracted IMF component and the residual component, respectively.
(3) Finally, all of the predicting results from the ANN, SVM and ANFIS models are combined to obtain the new output values, which are the final forecasting result for the groundwater level prediction.
Groundwater level time series Predicted values

The Hybrid EEMD-ANN, EEMD-SVM and EEMD-ANFIS Forecasting Models
The goal of this study was to improve the forecasting accuracy of groundwater levels by coupling the EEMD and data-driven models. A schematic diagram of the hybrid EEMD-ANN, EEMD-SVM and EEMD-ANFIS models is illustrated in Figure 4. As shown in Figure 4, the process of running the three hybrid models consists of three main steps:

The Hybrid EEMD-ANN, EEMD-SVM and EEMD-ANFIS Forecasting Models
The goal of this study was to improve the forecasting accuracy of groundwater levels by coupling the EEMD and data-driven models. A schematic diagram of the hybrid EEMD-ANN, EEMD-SVM and EEMD-ANFIS models is illustrated in Figure 4. As shown in Figure 4, the process of running the three hybrid models consists of three main steps:  (2) Secondly, the three data driven models ANN, SVM and ANFIS are developed to make the corresponding prediction using extracted IMF component and the residual component, respectively.
(3) Finally, all of the predicting results from the ANN, SVM and ANFIS models are combined to obtain the new output values, which are the final forecasting result for the groundwater level prediction.
(2) Secondly, the three data driven models ANN, SVM and ANFIS are developed to make the corresponding prediction using extracted IMF component and the residual component, respectively.
(3) Finally, all of the predicting results from the ANN, SVM and ANFIS models are combined to obtain the new output values, which are the final forecasting result for the groundwater level prediction.

Study Site and Data Preprocessing
The study area was located in Florida, United States. The annual mean temperature in the study area was 23.36 • C and the average annual mean precipitation was 980.61 mm in the past 16 years. All the observed data (including the monthly mean precipitation, monthly mean temperature, monthly maximum temperature, monthly minimum temperature, monthly mean lake level and monthly mean groundwater level) were collected to forecast the future groundwater fluctuations at well M1255 and STL185. For the obtained data from the past 16 years (from 1997 to 2012), data from the first 14 years was used for training and the last two years of data was used for validation. Meanwhile, the IMF components and one residual component of the groundwater level time-series were used to predict future groundwater fluctuations based on the EEMD. The results of the groundwater level prediction using the EEMD were analyzed and compared with the results obtained not using the EEMD. The locations of the observed wells in the study area are shown in Figure 5. Figure 6a-d illustrates the monthly mean precipitation, monthly mean temperature, monthly mean lake level, and monthly maximum temperature. Figure 7a,b illustrates the monthly groundwater levels at well M1255 and STL185.
Water 2018, 10, x FOR PEER REVIEW 7 of 20

Study Site and Data Preprocessing
The study area was located in Florida, United States. The annual mean temperature in the study area was 23.36 °C and the average annual mean precipitation was 980.61 mm in the past 16 years. All the observed data (including the monthly mean precipitation, monthly mean temperature, monthly maximum temperature, monthly minimum temperature, monthly mean lake level and monthly mean groundwater level) were collected to forecast the future groundwater fluctuations at well M1255 and STL185. For the obtained data from the past 16 years (from 1997 to 2012), data from the first 14 years was used for training and the last two years of data was used for validation. Meanwhile, the IMF components and one residual component of the groundwater level time-series were used to predict future groundwater fluctuations based on the EEMD. The results of the groundwater level prediction using the EEMD were analyzed and compared with the results obtained not using the EEMD. The locations of the observed wells in the study area are shown in Figure 5. Figure 6a-d illustrates the monthly mean precipitation, monthly mean temperature, monthly mean lake level, and monthly maximum temperature. Figure 7a,b illustrates the monthly groundwater levels at well M1255 and STL185.
To eliminate the dimensional differences, the time-series data was normalized using the following equation before the training stage.
where i x is the normalized data, x is the time-series data, min x is the minimum value and max x is the maximum value.    The partial autocorrelation function (PACF), as one of the statistical methods, was frequently used to select suitable models [10]. The PACF of well site M1255 from lag-0 to lag-18 is shown in Figure 8. Figure 8 suggests a significant correlation up to the lag-1 month for this time series within the confidence interval. The partial autocorrelation coefficients indicated that there was a one-month time lag for the monthly groundwater level time series as the input vector of the ANN, SVM and ANFIS models.  The partial autocorrelation function (PACF), as one of the statistical methods, was frequently used to select suitable models [10]. The PACF of well site M1255 from lag-0 to lag-18 is shown in Figure 8. Figure 8 suggests a significant correlation up to the lag-1 month for this time series within the confidence interval. The partial autocorrelation coefficients indicated that there was a one-month time lag for the monthly groundwater level time series as the input vector of the ANN, SVM and ANFIS models. To eliminate the dimensional differences, the time-series data was normalized using the following equation before the training stage.
where x i is the normalized data, x is the time-series data, x min is the minimum value and x max is the maximum value. The partial autocorrelation function (PACF), as one of the statistical methods, was frequently used to select suitable models [10]. The PACF of well site M1255 from lag-0 to lag-18 is shown in Figure 8. Figure 8 suggests a significant correlation up to the lag-1 month for this time series within the confidence interval. The partial autocorrelation coefficients indicated that there was a one-month time lag for the monthly groundwater level time series as the input vector of the ANN, SVM and ANFIS models.

Performance Criteria
Five statistical performance evaluation parameters-correlation coefficient (R), normalized mean square error (NMSE), root mean squared error (RMSE), Nash-Sutcliffe efficiency coefficient (NS) and Akaike information criteria (AIC)-were used to estimate the performance of the models proposed in this study.
Correlation coefficient (R): Normalized mean square error (NMSE): Root mean squared error (RMSE): Nash-Sutcliffe efficiency coefficient (NS): Akaike information criteria (AIC): where i O is the observed groundwater level value, i P is the predicted groundwater level value, O is the average of the observed groundwater level value, P is the average of the predicted groundwater level value, N is the number of input samples, and k is the number of free parameters used in those models. Suppose k is 0 in this paper. The Akaike information criteria (AIC) was employed to select the best time series mode [46]. The best fit between the predicted value and observed value would have R = 1, NMSE = 0, RMSE = 0, NS = 1, AIC = −∞ , respectively.

Decomposing, Modeling and Application
The time series of the original groundwater levels can be decomposed into several independent IMFs and one residue using the EEMD technique. The results are shown in Figures 9 and 10. It is

Performance Criteria
Five statistical performance evaluation parameters-correlation coefficient (R), normalized mean square error (NMSE), root mean squared error (RMSE), Nash-Sutcliffe efficiency coefficient (NS) and Akaike information criteria (AIC)-were used to estimate the performance of the models proposed in this study.
Correlation coefficient (R): Normalized mean square error (NMSE): Root mean squared error (RMSE): Nash-Sutcliffe efficiency coefficient (NS): Akaike information criteria (AIC): where O i is the observed groundwater level value, P i is the predicted groundwater level value, O is the average of the observed groundwater level value, P is the average of the predicted groundwater level value, N is the number of input samples, and k is the number of free parameters used in those models. Suppose k is 0 in this paper. The Akaike information criteria (AIC) was employed to select the best time series mode [46]. The best fit between the predicted value and observed value would have R = 1, NMSE = 0, RMSE = 0, NS = 1, AIC = −∞, respectively.

Decomposing, Modeling and Application
The time series of the original groundwater levels can be decomposed into several independent IMFs and one residue using the EEMD technique. The results are shown in Figures 9 and 10. It is noticeable that the two time series of the groundwater levels were resolved into six independent IMF and one residue component. Six independent IMFs were successively displayed from the highest to the lowest frequency, respectively. The extracted IMF components and the residual component as input variables were used to forecast the groundwater level. EEMD, as a new noise-assisted data analysis method, overcame the scale separation problem, avoided the mode mixing problem and conserved the physical uniqueness of the decomposition [19]. Therefore, the decomposition could be helpful to transform non-linear and non-stationary time series to independent IMF components from high to low frequency and could be useful to improve the forecast accuracy. noticeable that the two time series of the groundwater levels were resolved into six independent IMF and one residue component. Six independent IMFs were successively displayed from the highest to the lowest frequency, respectively. The extracted IMF components and the residual component as input variables were used to forecast the groundwater level. EEMD, as a new noise-assisted data analysis method, overcame the scale separation problem, avoided the mode mixing problem and conserved the physical uniqueness of the decomposition [19]. Therefore, the decomposition could be helpful to transform non-linear and non-stationary time series to independent IMF components from high to low frequency and could be useful to improve the forecast accuracy.  noticeable that the two time series of the groundwater levels were resolved into six independent IMF and one residue component. Six independent IMFs were successively displayed from the highest to the lowest frequency, respectively. The extracted IMF components and the residual component as input variables were used to forecast the groundwater level. EEMD, as a new noise-assisted data analysis method, overcame the scale separation problem, avoided the mode mixing problem and conserved the physical uniqueness of the decomposition [19]. Therefore, the decomposition could be helpful to transform non-linear and non-stationary time series to independent IMF components from high to low frequency and could be useful to improve the forecast accuracy.  The selection of input variables depended on the availability of the inputs and historically observed records. Several combinations of the exogenous factors were applied in the ANN, SVM and ANFIS models for groundwater level prediction. The input variables in the ANN model were the monthly time-series data of the precipitation, temperature (maximum, mean and minimum), and groundwater level. ANN+L representing the lake level was added to input variables in the ANN model. EEMD-ANN represents that the input variables were IMF components and one residual component. The input variables of the SVM and ANFIS models were similar to that of the ANN models (see Tables 1-3).

EEMD-ANN and ANN Models
For the investigation of the effects of the input structure and input data on the performance of the ANN model, three input structures were considered, as shown in Table 1. The results showed that the input variables of the EEMD-ANN at two sites has the minimum RMSE and AIC values in the validation. Meanwhile, the input variables of EEMD-ANN had the maximum R value in the validation. At site M1255, the minimum RMSE and AIC of EEMD-ANN were 0.329 and −53.350, respectively; at site STL185, the minimum RMSE and AIC of EEMD-ANN were 0.646 and −20.979, respectively. The maximum R of EEMD-ANN in the validation was 0.926 at site M1255 and 0.891 at site STL185. Table 1 suggested that R values of the four input variables of ANN+L at two sites were better than that of three input variables of ANN in the validation stage. The R of ANN+L at site M1255 was 0.754. The R value for ANN+L applied at site STL185 was 0.839. Figure 11 shows the observed and predicted groundwater levels using ANN, ANN+L and EEMD-ANN in the validation stage. The selection of input variables depended on the availability of the inputs and historically observed records. Several combinations of the exogenous factors were applied in the ANN, SVM and ANFIS models for groundwater level prediction. The input variables in the ANN model were the monthly time-series data of the precipitation, temperature (maximum, mean and minimum), and groundwater level. ANN+L representing the lake level was added to input variables in the ANN model. EEMD-ANN represents that the input variables were IMF components and one residual component. The input variables of the SVM and ANFIS models were similar to that of the ANN models (see Tables 1-3).

EEMD-ANN and ANN Models
For the investigation of the effects of the input structure and input data on the performance of the ANN model, three input structures were considered, as shown in Table 1. The results showed that the input variables of the EEMD-ANN at two sites has the minimum RMSE and AIC values in the validation. Meanwhile, the input variables of EEMD-ANN had the maximum R value in the validation. At site M1255, the minimum RMSE and AIC of EEMD-ANN were 0.329 and −53.350, respectively; at site STL185, the minimum RMSE and AIC of EEMD-ANN were 0.646 and −20.979, respectively. The maximum R of EEMD-ANN in the validation was 0.926 at site M1255 and 0.891 at site STL185. Table 1 suggested that R values of the four input variables of ANN+L at two sites were better than that of three input variables of ANN in the validation stage. The R of ANN+L at site M1255 was 0.754. The R value for ANN+L applied at site STL185 was 0.839. Figure 11 shows the observed and predicted groundwater levels using ANN, ANN+L and EEMD-ANN in the validation stage.

EEMD-SVM and SVM Models
The same input structures were introduced to the SVM model, as shown in Table 2. The results showed that the input variables of EEMD-SVM at site M1255 and the three input variables of SVM at site STL185 had the minimum RMSE and AIC values in the validation. The minimum RMSE and AIC

EEMD-SVM and SVM Models
The same input structures were introduced to the SVM model, as shown in Table 2. The results showed that the input variables of EEMD-SVM at site M1255 and the three input variables of SVM at site STL185 had the minimum RMSE and AIC values in the validation. The minimum RMSE and AIC of EEMD-SVM at site M1255 were 0.315 and 221255.373, respectively. The minimum RMSE and AIC of SVM at site STL185 were 0.755 and −13.512, respectively. Table 2 suggests that the maximum R value of EEMD-SVM was 0.885 at site M1255 and 0.860 at site STL185. The R values of the four input variables of SVM+L at the two sites were higher than or equal to those of the three input variables of SVM in the validation stage. The R values of SVM+L at site M1255 were 0.758 and 0.730 in the training and validation stage, respectively. The R values of SVM+L at site STL185 were 0.858 and 0.845 in the training and validation stage, respectively. Figure 12 shows the observed and predicted groundwater levels using SVM, SVM+L and EEMD-SVM in the validation stage.  Figure 12 shows the observed and predicted groundwater levels using SVM, SVM+L and EEMD-SVM in the validation stage.

EEMD-ANFIS and ANFIS Models
In the ANFIS model, the grid partition algorithm was used to generate a fuzzy inference system (FIS) structure from the training data [47]. The same input structures were introduced to the ANFIS model, the results are shown in Table 3. Table 3 suggested that the input variables of EEMD-ANFIS at two sites had the minimum RMSE and AIC values in the validation. The minimum RMSE and AIC of EEMD-ANFIS at site M1255 were 0.360 and −49.100, respectively. The minimum RMSE and AIC of EEMD-ANFIS at site STL185 were 0.606 and −24.035, respectively. The result showed that the maximum R value of EEMD-ANFIS in the validation was 0.926 at site M1255 and 0.909 at site STL185.
The R values of four input variables of ANFIS+L at two sites were greater than that of the three input variables of ANFIS in the validation stage. The R value of ANFIS+L at site M1255 was 0.799 in the validation. The R value of ANFIS+L at site STL185 was 0.910 in the validation. Figure 13 shows the observed and predicted groundwater levels using ANFIS, ANFIS+L and EEMD-ANFIS in the validation stage.

EEMD-ANFIS and ANFIS Models
In the ANFIS model, the grid partition algorithm was used to generate a fuzzy inference system (FIS) structure from the training data [47]. The same input structures were introduced to the ANFIS model, the results are shown in Table 3. Table 3 suggested that the input variables of EEMD-ANFIS at two sites had the minimum RMSE and AIC values in the validation. The minimum RMSE and AIC of EEMD-ANFIS at site M1255 were 0.360 and −49.100, respectively. The minimum RMSE and AIC of EEMD-ANFIS at site STL185 were 0.606 and −24.035, respectively. The result showed that the maximum R value of EEMD-ANFIS in the validation was 0.926 at site M1255 and 0.909 at site STL185. The R values of four input variables of ANFIS+L at two sites were greater than that of the three input variables of ANFIS in the validation stage. The R value of ANFIS+L at site M1255 was 0.799 in the validation. The R value of ANFIS+L at site STL185 was 0.910 in the validation. Figure 13 shows the observed and predicted groundwater levels using ANFIS, ANFIS+L and EEMD-ANFIS in the validation stage.

Comparison of EEMD-ANN, EEMD-SVM and EEMD-ANFIS
The best statistical parameter values of validation are shown in Tables 1-3, when the IMF components and one residual component were considered as the input variables. Also, based on the results of the statistical analysis, it was found that taking the lake level into account for the input variables led to better prediction in terms of the accuracy. Since the two well sites for data collection were close to the lake, groundwater-lake interaction could affect the groundwater level fluctuations. The results showed that the lake level should be considered as an input variable when exogenous factors were used to forecast the groundwater level for those well sites.
The analysis results from the model training and validation stages for two well sites are listed in Tables 1-3. At the training stage, the RMSE values in the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models for well M1255 were 0.234, 0.206 and 0.162 respectively; the RMSE values in these models for well STL185 were 0.281, 0.360 and 0.249 respectively. The RMSE value of the EEMD-ANFIS model was smaller than those of other two models in the training stage, which implied that the prediction ability of the EEMD-ANFIS model was better than that of the other two models for the given data. In the validation stage, the RMSE values for the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models for well M1255 were 0.329, 0.315 and 0.360 respectively; the RMSE of these models for well STL185 were 0.646, 1.009 and 0.606 respectively. For well M1255, the prediction result based on EEMD-ANFIS was close to that obtained from the other two models; the prediction result of EEMD-ANFIS was more accurate than that of the other two models for well STL185.
Generally, the modeling result is regarded as a perfect estimation when the NS criterion is equal to 1. If the NS criterion is higher than 0.8, the model can be recognized as effective and accurate [48]. The NS values for the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models in the training stage were all greater than 0.8. This indicated that the results from all these hybrid models are acceptable for forecasting groundwater levels. In the validation stage, the NS values for the EEMD-SVM model at site M1255 were greater than those for the EEMD-ANN and EEMD-ANFIS models. The NS values for the EEMD-ANFIS model at site STL185 were higher than those for the EEMD-ANN and EEMD-

Comparison of EEMD-ANN, EEMD-SVM and EEMD-ANFIS
The best statistical parameter values of validation are shown in Tables 1-3, when the IMF components and one residual component were considered as the input variables. Also, based on the results of the statistical analysis, it was found that taking the lake level into account for the input variables led to better prediction in terms of the accuracy. Since the two well sites for data collection were close to the lake, groundwater-lake interaction could affect the groundwater level fluctuations. The results showed that the lake level should be considered as an input variable when exogenous factors were used to forecast the groundwater level for those well sites.
The analysis results from the model training and validation stages for two well sites are listed in Tables 1-3. At the training stage, the RMSE values in the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models for well M1255 were 0.234, 0.206 and 0.162 respectively; the RMSE values in these models for well STL185 were 0.281, 0.360 and 0.249 respectively. The RMSE value of the EEMD-ANFIS model was smaller than those of other two models in the training stage, which implied that the prediction ability of the EEMD-ANFIS model was better than that of the other two models for the given data. In the validation stage, the RMSE values for the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models for well M1255 were 0.329, 0.315 and 0.360 respectively; the RMSE of these models for well STL185 were 0.646, 1.009 and 0.606 respectively. For well M1255, the prediction result based on EEMD-ANFIS was close to that obtained from the other two models; the prediction result of EEMD-ANFIS was more accurate than that of the other two models for well STL185.
Generally, the modeling result is regarded as a perfect estimation when the NS criterion is equal to 1. If the NS criterion is higher than 0.8, the model can be recognized as effective and accurate [48].
The NS values for the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models in the training stage were all greater than 0.8. This indicated that the results from all these hybrid models are acceptable for forecasting groundwater levels. In the validation stage, the NS values for the EEMD-SVM model at site M1255 were greater than those for the EEMD-ANN and EEMD-ANFIS models. The NS values for the EEMD-ANFIS model at site STL185 were higher than those for the EEMD-ANN and EEMD-SVM models in the validation stage. This indicated that EEMD-ANFIS and EEMD-SVM had an overall better estimation quality in comparison with the EEMD-ANN model. (See Tables 1-3).
When comparing the RMSE and R values of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models in the validation, the RMSE value of the EEMD-SVM model at site M1255 was less than that in both the EEMD-ANN and the EEMD-ANFIS model. The RMSE value of the EEMD-ANFIS model at site STL185 was higher than both the EEMD-ANN and the EEMD-SVM model. For the results for the two well sites, EEMD-ANFIS had higher R values compared to the other models. Obviously, the R and NS values of the EEMD-ANFIS and the EEMD-SVM model in the validation stage were greater than that of the EEMD-ANN model. Therefore, the EEMD-ANFIS and EEMD-SVM model could be good data-driven model in the validation stage. (See Tables 1-3). Figures 14 and 15 illustrate the coefficient of determination (R 2 ) values, corresponding to the predicted values in the scatter plots at the M1255 and STL185 observation wells, using the EEMD-ANN, EEMD-SVM, EEMD-ANFIS, ANN, SVM, ANFIS, ANN+L, SVM+L and ANFIS+L models. The scatter plots revealed the relationships between the predicted and observed groundwater levels for two observation wells. It can be seen clearly from the scatter plots that the EEMD-ANFIS model forecast the groundwater levels with less scatter for the two observed wells. Figures 14 and 15 show that EEMD-ANFIS had the best fit line compared to the other models. Figure 16a,b show the forecast groundwater levels versus observed groundwater levels using all of the models in the training stage and the validation stage.  Tables 1-3). When comparing the RMSE and R values of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models in the validation, the RMSE value of the EEMD-SVM model at site M1255 was less than that in both the EEMD-ANN and the EEMD-ANFIS model. The RMSE value of the EEMD-ANFIS model at site STL185 was higher than both the EEMD-ANN and the EEMD-SVM model. For the results for the two well sites, EEMD-ANFIS had higher R values compared to the other models. Obviously, the R and NS values of the EEMD-ANFIS and the EEMD-SVM model in the validation stage were greater than that of the EEMD-ANN model. Therefore, the EEMD-ANFIS and EEMD-SVM model could be good data-driven model in the validation stage. (See Tables 1-3). Figures 14 and 15 illustrate the coefficient of determination (R 2 ) values, corresponding to the predicted values in the scatter plots at the M1255 and STL185 observation wells, using the EEMD-ANN, EEMD-SVM, EEMD-ANFIS, ANN, SVM, ANFIS, ANN+L, SVM+L and ANFIS+L models. The scatter plots revealed the relationships between the predicted and observed groundwater levels for two observation wells. It can be seen clearly from the scatter plots that the EEMD-ANFIS model forecast the groundwater levels with less scatter for the two observed wells. Figures 14 and 15 show that EEMD-ANFIS had the best fit line compared to the other models. Figure 16a,b show the forecast groundwater levels versus observed groundwater levels using all of the models in the training stage and the validation stage.    According to the analysis, the R 2 value for the data at site M1255 indicated that EEMD-ANFIS performed better than the other models, although the RMSE value for the EEMD-SVM was less than that for EEMD-ANFIS. The R 2 value of the three hybrid models (i.e., EEMD coupled) were better than that of the other models. Therefore, both the EEMD-ANFIS and EEMD-SVM models can be considered good data-driven models at site M1255. Figure 15 shows that the R 2 value of EEMD-ANFIS at site STL185 was nearly equal to that of ANFIS+L and was a bit higher than that of EEMD-ANN. The RMSE value for EEMD-ANFIS was less than that for EEMD-ANN and ANFIS+L. Thus, the EEMD-ANFIS model can be considered the best estimation model at site STL185. For the two observation sites, the forecast results obtained from the three hybrid models were suggested to have better quality compared to those not coupled with EEMD.
Boxplot is importantly used to check whether the data-driven models are able to forecast these variations and corresponding prediction errors. It intuitively depicts the quartile values for the prediction error of groundwater data. Figures 17 and 18 display that the median value of the training errors is close to zero, which indicates the good performance of the data-driven models in the training stage in terms of the efficiency. Figures 17 and 18 show the comparison of errors between the results obtained by the EEMD-ANN, EEMD-SVM, EEMD-ANFIS, ANN, SVM, ANFIS, ANN+L, SVM+L and ANFIS+L models in the training period and the validation period. The results of the boxplots indicated that EEMD-ANFIS was the most accurate model in the training stage. In the validation stage, the result of the EEMD-ANFIS for well M1255 was close to that of the EEMD-SVM and EEMD-ANN; while in the prediction of the groundwater level for well STL185, the EEMD-ANFIS performed more accurately than the EEMD-SVM and EEMD-ANN.
Overall, the performance of the EEMD-ANFIS, EEMD-SVM and EEMD-ANN models was superior to the other models. The results of the ANFIS+L, SVM+L and ANN+L model were better 22   According to the analysis, the R 2 value for the data at site M1255 indicated that EEMD-ANFIS performed better than the other models, although the RMSE value for the EEMD-SVM was less than that for EEMD-ANFIS. The R 2 value of the three hybrid models (i.e., EEMD coupled) were better than that of the other models. Therefore, both the EEMD-ANFIS and EEMD-SVM models can be considered good data-driven models at site M1255. Figure 15 shows that the R 2 value of EEMD-ANFIS at site STL185 was nearly equal to that of ANFIS+L and was a bit higher than that of EEMD-ANN. The RMSE value for EEMD-ANFIS was less than that for EEMD-ANN and ANFIS+L. Thus, the EEMD-ANFIS model can be considered the best estimation model at site STL185. For the two observation sites, the forecast results obtained from the three hybrid models were suggested to have better quality compared to those not coupled with EEMD.
Boxplot is importantly used to check whether the data-driven models are able to forecast these variations and corresponding prediction errors. It intuitively depicts the quartile values for the prediction error of groundwater data. Figures 17 and 18 display that the median value of the training errors is close to zero, which indicates the good performance of the data-driven models in the training stage in terms of the efficiency. Figures 17 and 18 show the comparison of errors between the results obtained by the EEMD-ANN, EEMD-SVM, EEMD-ANFIS, ANN, SVM, ANFIS, ANN+L, SVM+L and ANFIS+L models in the training period and the validation period. The results of the boxplots indicated that EEMD-ANFIS was the most accurate model in the training stage. In the validation stage, the result of the EEMD-ANFIS for well M1255 was close to that of the EEMD-SVM and EEMD-ANN; while in the prediction of the groundwater level for well STL185, the EEMD-ANFIS performed more accurately than the EEMD-SVM and EEMD-ANN.
Overall, the performance of the EEMD-ANFIS, EEMD-SVM and EEMD-ANN models was superior to the other models. The results of the ANFIS+L, SVM+L and ANN+L model were better than that of ANFIS, SVM and ANN in terms of R. It indicated that lake level fluctuations as an input variable is important in the prediction of the groundwater level in the near-lake area. Comparing the results of the EEMD-ANFIS, EEMD-ANN and EEMD-SVM models, the three models had their own advantages and disadvantages. The selection of the prediction model should balance the effects and benefits of the statistical parameters (e.g., R and RMSE) in both the training period and the validation period. The prediction results of the EEMD-ANFIS were close to that of the EEMD-SVM and EEMD-ANN models at site M1255; the prediction results of the EEMD-ANFIS were a bit more accurate than that of the EEMD-SVM and EEMD-ANN models at site STL185. The results of the application suggested that the EEMD-ANFIS, EEMD-SVM and EEMD-ANN models were feasible and effective. Also, EEMD can be used to improve the accuracy of predicting nonlinear and nonstationary time series. The results in this study are consistent with those acquired by [22,26,29].
Water 2018, 10, x FOR PEER REVIEW 16 of 20 than that of ANFIS, SVM and ANN in terms of R. It indicated that lake level fluctuations as an input variable is important in the prediction of the groundwater level in the near-lake area. Comparing the results of the EEMD-ANFIS, EEMD-ANN and EEMD-SVM models, the three models had their own advantages and disadvantages. The selection of the prediction model should balance the effects and benefits of the statistical parameters (e.g., R and RMSE) in both the training period and the validation period. The prediction results of the EEMD-ANFIS were close to that of the EEMD-SVM and EEMD-ANN models at site M1255; the prediction results of the EEMD-ANFIS were a bit more accurate than that of the EEMD-SVM and EEMD-ANN models at site STL185. The results of the application suggested that the EEMD-ANFIS, EEMD-SVM and EEMD-ANN models were feasible and effective. Also, EEMD can be used to improve the accuracy of predicting nonlinear and nonstationary time series. The results in this study are consistent with those acquired by [22,26,29].

Conclusions
The reliable and accurate estimation of groundwater level fluctuation is essential in order to manage water resources and improve water-use efficiency. In this study, the prediction capability of the ANN, SVM and ANFIS models based on ensemble empirical mode decomposition (EEMD) were investigated using monthly groundwater level data collected at the M1255 and STL185 observation wells. The statistical parameters R, NMSE, RMSE, NS and AIC were used to assess the performance of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models. The results from the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models were analyzed and compared with the results from the ANN, SVM and ANFIS models. The values of the statistical parameters indicated that the EEMD-ANN, EEMD-SVM, and EEMD-ANFIS models achieved better prediction results than the ANN, SVM and ANFIS. The average R value of three hybrid models was higher than that of the ANN, SVM and ANFIS models, and the average RMSE value of these hybrid models was less than that of the ANN, SVM and ANFIS models.
The results in this study suggested that EEMD can effectively enhance predicting accuracy. The proposed EEMD could significantly improve the performance of the ANN, SVM and ANFIS date-driven models in groundwater level forecasting. The proposed three hybrid models based on EEMD had several obvious advantages: (a) it was convenient and effective to combine the EEMD with the ANN, SVM and ANFIS to forecast the nonstationary and nonlinear groundwater level fluctuations; (b) time series data on the groundwater level was only required in the hybrid models and exogenous factors affecting the groundwater level do not need to be considered in the research area; and (c) the prediction results of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models were more accurate when using the groundwater level time series decomposition. Therefore, this study supported the validity and applicability of the EEMD-ANN, EEMD-SVM and EEMD-ANFIS models in the prediction of groundwater levels. The results from this research would be beneficial for sustainable water resource management.
Author Contributions: Y.G. wrote the manuscript; Z.W., G.X. and Z.Z. revised and proofread the manuscript.