Application of a Novel Hybrid Wavelet-ANFIS/Fuzzy C-Means Clustering Model to Predict Groundwater Fluctuations

: In order to optimize the management of groundwater resources, accurate estimates of groundwater level (GWL) ﬂuctuations are required. In recent years, the use of artiﬁcial intelligence methods based on data mining theory has increasingly attracted attention. The goal of this research is to evaluate and compare the performance of adaptive network-based fuzzy inference system (ANFIS) and Wavelet-ANFIS models based on FCM for simulation/prediction of monthly GWL in the Maragheh plain in northwestern Iran. A 22-year dataset (1996–2018) including hydrological parameters such as monthly precipitation (P) and GWL from 25 observation wells was used as models input data. To improve the prediction accuracy of hybrid Wavelet-ANFIS model, different mother wavelets and different numbers of clusters and decomposition levels were investigated. The new hybrid model with Sym4-mother wavelet, two clusters and a decomposition level equal to 3 showed the best performance. The maximum values of R2 in the training and testing phases were 0.997 and 0.994, respectively, and the best RMSE values were 0.05 and 0.08 m, respectively. By comparing the results of the ANFIS and hybrid Wavelet-ANFIS models, it can be deduced that a hybrid model is an acceptable method in modeling of GWL because it employs both the wavelet transform and FCM clustering technique.


Introduction
Hydrological changes, also linked to global climate change, will have devastating effects on surface water and groundwater resources in many areas of the world [1]. Many areas are already facing water shortages such as the Middle East, including Iran. Iran has arid and semi-arid climates with limited surface water resources, so the main water resource is groundwater in most regions [2,3]. Recent droughts, coupled with irregular groundwater withdrawal through deep and semi-deep wells, have caused a severe drop in groundwater levels (GWL) and consequently water quality degradation, soil erosion, drying of wells, reducing river discharge, increasing cost of pumping, and desertification [4]. In this regard, an effective water resources management is vital. The monitoring, simulation, and prediction of GWL help water planners apply appropriate measures for water resources management. This can be done by a field measurement of GWL using borehole data which is, however, expensive and time consuming. In recent years, the use of indirect methods such as artificial intelligence (AI) and machine learning (ML) models has been developed to analyze the trends of hydrological parameters. The latter have been used in many hydro(geo)logical studies, specifically in GWL predictions [5][6][7][8][9][10][11]. Finding suitable models

Data Description
In order to simulate and predict the GWL fluctuations in the Maragheh plain, the monthly averages of rainfall and GWL data were collected from the State Government. The observed precipitation and GWL (Table 1) are available for the years 1996 to 2018 (261 months). The locations of these wells with their Thiessen polygons are shown in Figure 2.

Data Description
In order to simulate and predict the GWL fluctuations in the Maragheh plain, the monthly averages of rainfall and GWL data were collected from the State Government. The observed precipitation and GWL (Table 1) are available for the years 1996 to 2018 (261 months). The locations of these wells with their Thiessen polygons are shown in Figure 2.

Completion of Groundwater Level Missing Values
Among 25 observation wells, only 18 wells had complete time series of raw GWL-data during the period 1996 to 2018. The missing values of monthly GWL for seven wells (see Table 2) were estimated by inverse distance weighting (IDW) interpolation. The credibility of the estimated values was checked by comparison with previous month values. Subsequently, the GWL-data were average-weighted using the (25) Thiessen polygon zones to create weight-hydrographs of the GWL-fluctuations ( Figure 3).

Completion of Groundwater Level Missing Values
Among 25 observation wells, only 18 wells had complete time series of raw GWLdata during the period 1996 to 2018. The missing values of monthly GWL for seven wells (see Table 2) were estimated by inverse distance weighting (IDW) interpolation. The credibility of the estimated values was checked by comparison with previous month values. Subsequently, the GWL-data were average-weighted using the (25) Thiessen polygon zones to create weight-hydrographs of the GWL-fluctuations ( Figure 3).

Training and Testing of GWL Prediction Models
The ANFIS and Wavelet-ANFIS models were trained based on monthly-observed GWL as well as precipitation (P) data time series. The following input-output prediction model for GWL at time t as a function of GWLt−1 and Pt−j was applied:

Training and Testing of GWL Prediction Models
The ANFIS and Wavelet-ANFIS models were trained based on monthly-observed GWL as well as precipitation (P) data time series. The following input-output prediction model for GWL at time t as a function of GWL t−1 and P t−j was applied: i.e., it is assumed that the groundwater levels in month t (GWLt) are correlated with up to t-i previous months 'groundwater levels and precipitation P t−j . Although Equation (1) is basically just a simple transfer model, it has some physical basis, as groundwater data usually exhibit some persistence, as illustrated in Figure 3, given that groundwater is recharged by rainfall with a delay [40,41]. Autocorrelation analysis of the GWL-time series and cross correlation analysis of the GWL-with the precipitation time series P were performed, respectively. The continuous autocorrelation function Rxx of a time series x is defined as [41]: where T is the period of observation and l is the lag time. For discrete data, Equation (2) is replaced by its discrete homologue which for x = GWL, reads where n (=264) is the number of observed months. The cross-correlation function Rxy between two time series x and y is defined similarly to Equations (2) and (3), so that the discrete cross correlation between GWL and P at lag l reads as: Figure 4 shows that GWL t is mostly correlated with GWL t−1,t−2,t−3 and P t−1,t−2,t−3,t−4 . Therefore, the final input-output Wavelet-ANFIS and ANFIS and SVM model Equation (1) can be written as: Atmosphere 2021, 12, x FOR PEER REVIEW 5 of 15 i.e., it is assumed that the groundwater levels in month t (GWLt) are correlated with up to t-i previous months 'groundwater levels and precipitation Pt−j. Although Equation (1) is basically just a simple transfer model, it has some physical basis, as groundwater data usually exhibit some persistence, as illustrated in Figure 3, given that groundwater is recharged by rainfall with a delay [40,41]. Autocorrelation analysis of the GWL-time series and cross correlation analysis of the GWL-with the precipitation time series P were performed, respectively. The continuous autocorrelation function Rxx of a time series x is defined as [41]: where T is the period of observation and l is the lag time. For discrete data, Equation (2) is replaced by its discrete homologue which for x = GWL, reads where n (=264) is the number of observed months. The cross-correlation function Rxy between two time series x and y is defined similarly to Equations (2) and (3), so that the discrete cross correlation between GWL and P at lag l reads as: Figure 4 shows that GWLt is mostly correlated with GWLt−1,t−2,t−3 and Pt−1,t−2,t−3,t−4. Therefore, the final input-output Wavelet-ANFIS and ANFIS and SVM model Equation (1) can be written as: For the training and testing of the models the cross-validation method has been used, i.e., the 261-values observed average GWL-and P-data series are divided into a training set of 183 months and testing set of 78 months for prediction with the trained model. Thus, the length of the training data set is 70% and the test (prediction) data makes up the remaining 30%.

Adaptive Neuro Fuzzy Inference System (ANFIS)
The adaptive neuro fuzzy inference system (ANFIS), first introduced by [13], is a combination of an adaptive neural network and a fuzzy inference system. This model has also an adequate in training, producing and categorizing and has several advantages over classical ANN, such as the capability of handling large amounts of data, dynamic and nonlinear systems modeling, easy use, and efficient computing, while still exhibiting high estimation and prediction accuracies [42]. Combining at the same time the benefits and capabilities of neural network architecture methods and fuzzy logic, ANFIS uses a hybrid For the training and testing of the models the cross-validation method has been used, i.e., the 261-values observed average GWL-and P-data series are divided into a training set of 183 months and testing set of 78 months for prediction with the trained model. Thus, the length of the training data set is 70% and the test (prediction) data makes up the remaining 30%.

Adaptive Neuro Fuzzy Inference System (ANFIS)
The adaptive neuro fuzzy inference system (ANFIS), first introduced by [13], is a combination of an adaptive neural network and a fuzzy inference system. This model has also an adequate in training, producing and categorizing and has several advantages over classical ANN, such as the capability of handling large amounts of data, dynamic and nonlinear systems modeling, easy use, and efficient computing, while still exhibiting high estimation and prediction accuracies [42]. Combining at the same time the benefits and capabilities of neural network architecture methods and fuzzy logic, ANFIS uses a hybrid approach of the classical gradient descent procedure and systematic back-propagation and thus avoids the "trap" of the error function in a local minimum, typical of the former optimization method [43]. The adaptive network based on the Sugeno fuzzy inference model provides a deterministic system of output equations and is thus a useful approach for parameter estimation [44]. The ANFIS model consists of five layers and is given in Figure 5 [13]. The first layer is called input layer; O 1,i is the output of the ith node of the layer 1.
where X 1,2 is the input node i and A i (or B i ) is a linguistic label associated with this node. Therefore O 1,i is the membership grade of a fuzzy set (A 1 , A 2 , B 1 , B 2 ). The second layer consists of rule nodes with AND and/or OR operators. The output (O 2,i ) is the product of all the incoming signals.
Atmosphere 2021, 12, x FOR PEER REVIEW 6 of 15 approach of the classical gradient descent procedure and systematic back-propagation and thus avoids the "trap" of the error function in a local minimum, typical of the former optimization method [43]. The adaptive network based on the Sugeno fuzzy inference model provides a deterministic system of output equations and is thus a useful approach for parameter estimation [44]. The ANFIS model consists of five layers and is given in Figure 5 [13]. The first layer is called input layer; O1,i is the output of the ith node of the layer 1.
where X1,2 is the input node i and Ai (or Bi) is a linguistic label associated with this node. Therefore O1,i is the membership grade of a fuzzy set (A1, A2, B1, B2). The second layer consists of rule nodes with AND and/or OR operators. The output (O2,i) is the product of all the incoming signals. Figure 5. Architecture of ANFIS [11].
The third layer outputs (O3,i) are called normalized firing strengths.
The fourth layer, called consequent node, is a standard perceptron. Every node i in this layer is an adaptive node with a node function: where (pi, qi, ri) is the parameter set of this node; These are referred to as the consequent parameters. Finally, the fifth layer is called the output layer, which computes the overall output as the summation of all incoming signals. , = overall output = = ∑ (10) Figure 5. Architecture of ANFIS [11].
The third layer outputs (O 3,i ) are called normalized firing strengths.
The fourth layer, called consequent node, is a standard perceptron. Every node i in this layer is an adaptive node with a node function: Atmosphere 2021, 12, 9 7 of 15 where (p i , q i , r i ) is the parameter set of this node; These are referred to as the consequent parameters. Finally, the fifth layer is called the output layer, which computes the overall output as the summation of all incoming signals.
2.6. Fuzzy C-Means (FCM) Clustering Method ANFIS uses neural network training methods and employs fuzzy logic to fit a relationship between some inputs space to output space [13]. The problem with the former is that as the number of inputs increases, the number of rules rises rapidly. In other words, the fuzzy rules increase exponentially [22]. In order to avoid the dimensionality problem in generating too large rule-bases in ANFIS, clustering algorithms are used in FIS models to generate fewer fuzzy rules. It should be noted that clustering methods provide information to obtain initial rules and fuzzy structure for generating the fuzzy inference system [45]. Then, the fuzzy inference system based on these clustering methods can provide a powerful model for establishing relationships between inputs and the output space. The fuzzy C-means (FCM) clustering method was initially introduced by Bezdek in 1973 [46,47] and this approach is known as an improvement and modification of the K-means clustering. In this approach, each data point belongs partly to all clusters, however, with different membership degrees that can vary between 0 and 1. The objective function in the FCM is minimizing the total intra-cluster variance, i.e., the summed square error function within clusters (j m ) [48]: where U = [u ij ] c×n is the fuzzy partition matrix of c clusters and n data, V = (v 1 , v 2 , . . . , v i ) are the centroids of the clusters, u ij is the partial membership degree of data x j in cluster . . , n), m is the weighting exponent on each fuzzy membership (which is equal to 2 in this study), x j = (x 1 , x 2 , . . . , x n ) is the dataset, v i is the initial value for the cluster center, and d(x j , v i ) is the Euclidean distance between the x j data and the cluster center of the i th cluster v i , i.e., x j −v i . Equation (11) describes a non-linear optimization problem which can be solved by iterative minimization. The centroid of each cluster is calculated by the partial derivative of Equation (11) with respect to V, then the partial membership degree of data is updated in each iteration by differentiation of the above equation with respect to U: The iteration process is stopped when the variance of intra-clusters in five iterations does not change more than the determined minimum improvement (in this study 10 −8 ), i.e., a minimum improvement in the FCM objective function (convergence) does not occur in five iterations and/or a maximum number of iterations (1000 iterations in this study) is exceeded. In the present study, the FCM algorithm was incorporated into (i.e., coded in) the ANFIS model. The procedure is sketched in Figure 6. The optimal value of the number of clusters was obtained based on trial and error. The RMSE values were calculated for different numbers in the training and test phases ( Figure 7). As shown in Figure 7, using two clusters resulted in the least error.

Wavelet Transformation
The Wavelet Transformation (WT) method is a mathematical tool in Fourier Transform, but seems a more effective tool than the Fourier transform (FT) for studying nonstationary signals [49]. WT is popular because of its multi resolution in time and frequency domain related to temporal data and has been developed to determine information that cannot be easily achieved from raw signals themselves. WT has the ability to consider different aspects of a time series such as trends, discontinuities, and breakpoints [50,51]. In the past, WT was widely used for many hydrological processes including GWL fluctuations, river flow forecast, fault classification or suspended sediments. The continuous wavelet transform (CWT) W(τ,s) of a signal x(t) with respect to a mother wavelet ψ(t) is defined as follows [52]: where τ is the translation parameter, s is the wavelet scale parameter, t is time, and * refers to the complex conjugate. However, the CWT is not often used for forecasting because it requires calculating wavelet coefficients at each scale which means more calculations and, in turn, more compute time. Instead, the discrete wavelet transform (DWT) requires less computation time and is simpler to develop [53]. Discrete wavelets have the following general form [49]:

Wavelet Transformation
The Wavelet Transformation (WT) method is a mathematical tool in Fourier Transform, but seems a more effective tool than the Fourier transform (FT) for studying nonstationary signals [49]. WT is popular because of its multi resolution in time and frequency domain related to temporal data and has been developed to determine information that cannot be easily achieved from raw signals themselves. WT has the ability to consider different aspects of a time series such as trends, discontinuities, and breakpoints [50,51]. In the past, WT was widely used for many hydrological processes including GWL fluctuations, river flow forecast, fault classification or suspended sediments. The continuous wavelet transform (CWT) W(τ,s) of a signal x(t) with respect to a mother wavelet ψ(t) is defined as follows [52]: where τ is the translation parameter, s is the wavelet scale parameter, t is time, and * refers to the complex conjugate. However, the CWT is not often used for forecasting because it requires calculating wavelet coefficients at each scale which means more calculations and, in turn, more compute time. Instead, the discrete wavelet transform (DWT) requires less computation time and is simpler to develop [53]. Discrete wavelets have the following general form [49]:

Wavelet Transformation
The Wavelet Transformation (WT) method is a mathematical tool in Fourier Transform, but seems a more effective tool than the Fourier transform (FT) for studying nonstationary signals [49]. WT is popular because of its multi resolution in time and frequency domain related to temporal data and has been developed to determine information that cannot be easily achieved from raw signals themselves. WT has the ability to consider different aspects of a time series such as trends, discontinuities, and breakpoints [50,51]. In the past, WT was widely used for many hydrological processes including GWL fluctuations, river flow forecast, fault classification or suspended sediments. The continuous wavelet transform (CWT) W(τ,s) of a signal x(t) with respect to a mother wavelet ψ(t) is defined as follows [52]: where τ is the translation parameter, s is the wavelet scale parameter, t is time, and * refers to the complex conjugate. However, the CWT is not often used for forecasting because it requires calculating wavelet coefficients at each scale which means more calculations and, in turn, more compute time. Instead, the discrete wavelet transform (DWT) requires less computation time and is simpler to develop [53]. Discrete wavelets have the following general form [49]: where m and n are integers that control the scale and time respectively, τ 0 is the location parameter that must be greater than 0, and S 0 is a specified fixed dilation step greater than 1. The DWT performs two functions viewed as high-pass and low-pass filters, through which the original time series are passed and then the original time series data are decomposed and divided into two parts, namely "approximation" and "details" [54]. These components explain behavior better and reveal more information about the process than the original time series. Therefore, they can help forecasting models predict with greater accuracy [55,56]. There are many wavelets that can be used as mother wavelets, For the choice of the discrete father-scaling/mother-wavelet filter functions in the DWT/MRA, popular ones are the Haar, the Daubechies wavelet db4 and the irregular wavelet Symlet (sym4), and these are those have been used here [38,57].

Hybrid Wavelet-ANFIS Model
As shown by [58], a combination of DWT and ANFIS can lead to greater accuracy of GWL prediction. In this regard, after the determination of important parameters, the ANFIS model with the FCM rule generator and the DWT-MRA-decomposed time series of the average GWL and P are integrated to form the Wavelet-ANFIS hybrid model. More specifically, several combinations of different time-lagged GWL-and P-data have been decomposed by haar, Db4, and sym4 data to 1, 2, 3, and 4 levels.

Results and Discussion
The performance of the ANFIS and Wavelet-ANFIS models were evaluated by calculating R 2 and RMSE statistical parameters. Shown in Table 3, as expected, the performance of the ANFIS model is lower than that of the Wavelet-ANFIS hybrid model, because the Wavelet-ANFIS hybrid model takes advantage of both ANFIS and WT, simultaneously. In particular, WT, or specifically, the DWTMRA, has the ability to consider various aspects of a time series such as trends, discontinuities and breakpoints, and this time-scale localization feature was used to provide decomposition of the input time series up to the 4 levels in the present application, allowing for a separation of the approximation and for providing details of the non-stationary noisy time series, thus enhancing the ANFIS, which uses a combination of fuzzification of the input through membership functions with the network-based algorithm and the use of the hybrid (back propagation tries and gradient descent) optimization method. Table 3 shows the minimum RMSE and the maximum R 2 in the training phase with 0.05 and 0.997 and in the test phase as 0.08 and 0.974, respectively, which are obtained with the hybrid wavelet ANFIS model and using the Sym4-mother wavelet (second input combination, Level Decomposition = 3). The best results of the Haar mother wavelets in the test phase has RMSE and R 2 values of 0.13 and 0.929, respectively. The db4 mother wavelets also exhibits the best results in the test phase, with corresponding values of 0.09 and 0.968. Comparison of the results of the db4 and sym4 mother wavelets revealed that, although there was no significant difference between the results of these two mother wavelets, the sym4 mother wavelets outperforms db4 mother wavelets slightly. The less accurate results are related to Haar mother wavelets since these are only orthogonal wavelet with a linear phase that cannot provide nonlinear shifts between the original and decomposed signals. By adding the precipitation parameter to the input data, the output accuracy increased in most simulations (on average, 10.2%), although it was more significant in the Wavelet-ANFIS models. In this research, the effect of the decomposition level on the accuracy of the output results of the Wavelet-ANFIS models was investigated. Among the different levels, the decomposition at level 3 had the highest accuracy. However, the accuracy of decomposition level 2 was also acceptable (with only a small difference compared to level 3).  Figure 8 shows the regression plots of the Wavelet-ANFIS-simulated over the observed GWL-data for the training and testing phases and in agreement with the high R 2 , the points are nearly perfectly lying on the "slope = 1" regression line. The GWL-time series observed and simulated with the best combination of the Wavelet-ANFIS model are presented in Figure 9. There is an acceptable correlation between the observed and simulated GWL data, and the effect of the applied delays on the input data, especially the piezometric head data, is clearly visible. All data-driven prediction methods are based on the idea that the random errors are drawn from a normal distribution and the hybrid Wavelet-ANFIS model is not an exception. Figure 10 illustrates that a-posteriori computed errors of the GWLs predicted with this model, follow indeed a normal distribution.
Exceptionally good and reliable predictions of the average GWL of the Maragheh plain using the proposed hybrid wavelet-ANFIS model are evident in all the above presented results. Although the average of GWL is the first and appropriate criterion for the evaluation of the GWL changes all over the plain, in order to understand the spatial variations across the piezometer network in the study plain, the GWL fluctuations in the 25 observation wells were simulated by the hybrid model. The calculated values of R 2 and RMSE for both the training and test phase are given in Figure 11. According to most R 2 values above 0.95 and most RMSE values less than 0.20 m, especially in the test phase, it can be argued that the model has an acceptable performance in almost all wells.
series observed and simulated with the best combination of the Wavelet-ANFIS model are presented in Figure 9. There is an acceptable correlation between the observed and simulated GWL data, and the effect of the applied delays on the input data, especially the piezometric head data, is clearly visible. All data-driven prediction methods are based on the idea that the random errors are drawn from a normal distribution and the hybrid Wavelet-ANFIS model is not an exception. Figure 10 illustrates that a-posteriori computed errors of the GWLs predicted with this model, follow indeed a normal distribution.     Exceptionally good and reliable predictions of the average GWL of the Maragheh plain using the proposed hybrid wavelet-ANFIS model are evident in all the above presented results. Although the average of GWL is the first and appropriate criterion for the evaluation of the GWL changes all over the plain, in order to understand the spatial vari- The results clearly indicate that simulating and predicting a vital parameter of the groundwater budget, namely GWL at time t, can be done accurately using the proposed hybrid model with values of GWL from the previous month and precipitation data as inputs. Although, this is a data driven model, choosing relevant inputs is important because, for example as illustrated in this study, the model considers the groundwater level changes that have a direct relationship with changes in groundwater storage (DS). Furthermore, the difference between volume of precipitation and DS is equal to the sum of other components of the groundwater budget equation, including irrigation water, domestic and industrial demands, as well as natural hydrological processes.

Conclusions
This study developed and evaluated hybrid Wavelet-ANFIS models based on the FCM method for groundwater level-(GWL) simulation and prediction, using the Maragheh plain, in northwestern Iran as a case study. In the developed models, GWL were predicted using antecedent GWL-data and precipitation data. A monthly GWL-data series recorded at 25 observation wells during a 22-year period was taken into account. However, given that the GWL data series were not complete, statistical data deficiencies were interpolated using the IDW method. After completing data pre-processing, various statistical analyses were carried out to investigate the time lag correlation of the GWL-time series it was established that that for yielding acceptable simulations of GWL values for a particular month, GWL-and P-data values from 3-4 months prior would be necessary. The 7-time lagged data were used directly as inputs to the ANFIS model to predict GWL, while in the hybrid Wavelet-ANFIS model, the GWL-and P-time series were decomposed to level 4 by means of DWT-MRA.
The performance of the developed hybrid model was examined at four decomposition levels. Input data were divided into two parts: one for simulation (training) and the other for prediction (testing) phases. Thus, in this study, the length of the training data set is 70% and the test (prediction) data makes up the remaining 30%. The best values of Figure 11. Results of GWL-simulation in all 25 wells (see Figure 2).
The results clearly indicate that simulating and predicting a vital parameter of the groundwater budget, namely GWL at time t, can be done accurately using the proposed hybrid model with values of GWL from the previous month and precipitation data as inputs. Although, this is a data driven model, choosing relevant inputs is important because, for example as illustrated in this study, the model considers the groundwater level changes that have a direct relationship with changes in groundwater storage (DS). Furthermore, the difference between volume of precipitation and DS is equal to the sum of other components of the groundwater budget equation, including irrigation water, domestic and industrial demands, as well as natural hydrological processes.

Conclusions
This study developed and evaluated hybrid Wavelet-ANFIS models based on the FCM method for groundwater level-(GWL) simulation and prediction, using the Maragheh plain, in northwestern Iran as a case study. In the developed models, GWL were predicted using antecedent GWL-data and precipitation data. A monthly GWL-data series recorded at 25 observation wells during a 22-year period was taken into account. However, given that the GWL data series were not complete, statistical data deficiencies were interpolated using the IDW method. After completing data pre-processing, various statistical analyses were carried out to investigate the time lag correlation of the GWL-time series it was established that that for yielding acceptable simulations of GWL values for a particular month, GWL-and P-data values from 3-4 months prior would be necessary. The 7-time lagged data were used directly as inputs to the ANFIS model to predict GWL, while in the hybrid Wavelet-ANFIS model, the GWL-and P-time series were decomposed to level 4 by means of DWT-MRA.
The performance of the developed hybrid model was examined at four decomposition levels. Input data were divided into two parts: one for simulation (training) and the other for prediction (testing) phases. Thus, in this study, the length of the training data set is 70% and the test (prediction) data makes up the remaining 30%. The best values of RMSE and R 2 were obtained 0.05 m and 0.997, respectively, in the training phase and of 0.08 m and 0.974 the testing (prediction) phase, which indicates that the hybrid Wavelet-ANFIS using the Symlet mother wavelet and decomposition level 3 performs best. According to the results of the hybrid Wavelet-ANFIS model, the model performance was acceptable in estimating GWL fluctuations all over the plain, such that the difference between the observed and simulated values was negligible, meaning the incorporation of the wavelet transforms in ANFIS increased the performance of ANFIS model, especially in the prediction phase. Another noteworthy characteristic of this coupled model in the simulated and predicted values of GWL fluctuations is the use of the FCM clustering method which generates fewer fuzzy rules, thereby overcoming the well-known data dimensionality problem.
Finally, it can be argued that the Wavelet-ANFIS-FCM-clustering model in terms of performance, especially, when, due to a lack of data, all physical processes such as surface-groundwater interactions, are not completely understood, presents an attractive method and reliable tool for groundwater level prediction under different water resources management scenarios. Nonetheless, there is a need to further evaluate the developed model in other regions or in the same plain but with different GWL-and P-time series. . The ground water data in Iran are not publicly available but can be requested from respective water authorities.

Conflicts of Interest:
The authors declare no conflict of interest.