Meticulous prediction of hydrological processes, especially water budget, has an individual importance in environmental management plans. On the other hand, conservation of groundwater, a fundamental resource in arid and semi-arid areas, needs to be considered as a great priority in development plans. Prediction of a groundwater budget utilizing artificial intelligence was the scope of this study. For this aim, the Azarshahr Plain aquifer, East Azerbaijan, Iran, was selected because of its great dependence on groundwater and the necessity of cognizance of its budget in future programs. The long-term fluctuations of the water table in 13 piezometers were simulated by a wavelet-based artificial neural network (WANN) hybrid model, and their statistical gaps were covered. Then, the modelled water table was predicted for the next 12 months using genetic programming. The results of simulation and prediction were assessed by performance evaluation criteria such as R^{2}, root mean squared error, mean absolute error and Nash–Sutcliffe efficiency. Thiessen polygons were then utilized, plotting the predicted unit hydrograph of the study area. The predicted water table from September 2012 to August 2013 revealed about 0.12 m depletion. Regarding the area of the Azarshahr Plain aquifer and its average storage coefficient, the aquifer budget will be reduced by about 0.3557 million cubic metres during this period.

## INTRODUCTION

Awareness of coming natural events, particularly hydrological processes, is a great challenge for environmental custodians, especially hydrologists. Notwithstanding their highly stochastic nature, the development of models capable of describing such complex phenomena is a growing area of research.

Groundwater as a major source of water supply for domestic, agricultural and industrial users and, of course, the main part of the hydrology cycle has a vital role in arid and semi-arid areas. In several such areas, much more groundwater is withdrawn than the recharge rate, leading to damaging environmental effects such as water level depletion, drying up of wells, abatement of water quality, amplified pumping prices and reduced well yields (Adamowski & Chan 2011). Effectively managing groundwater, predicting groundwater level ﬂuctuations, and quantifying these changes are strategic hydrological issues.

Recent years have seen a signiﬁcant rise in the number of scientiﬁc approaches applied to hydrologic modelling and forecasting, including the main popular ‘data-based’ or ‘data-driven’ approaches. Such modelling methods involve mathematical equations drawn not from the physical process in the watershed but from an analysis of simultaneous input and output time series (Solomatine & Ostfeld 2008). Meanwhile, it is becoming increasingly difficult to ignore the role of artificial intelligence (AI) in hydrological processes' prediction due to its efficiency in modeling complex physical processes based on certain data/information governing the process. There are numerous researches about groundwater level changes utilizing AI (e.g., Adamowski & Chan 2011; Kisi & Shiri 2012; Fallah-Mehdipour *et al.* 2013; Maheswaran & Khosa 2013; Moosavi *et al.* 2013; Nourani *et al.* 2015; Seo *et al.* 2015; Sivapragasam *et al.* 2015).

Nevertheless, AI is progressively being preferred primarily and many studies have already been reported in forecasting groundwater level changes; however, there is no direct forecasting of groundwater budget using this method. The low number of studies on groundwater budget modelling via AI demonstrates the need to consider groundwater and relevant issues.

Among various conceptual and black box models developed over the mentioned period, hybrid AI-based models have been among the most promising in simulating hydrologic processes (Nourani *et al.* 2014), and wavelet-AI is an example of these methods. On the other hand, genetic programming (GP) is an AI method that is based on the random iterative searching process to achieve an appropriate relationship between input and output. Conjugating the wavelet and GP methods can give an incredibly precise result in hydrological processes' prediction; for instance, Nourani *et al.* (2012) investigated the linkage of wavelet analysis to GP in constructing a hybrid model to detect seasonality patterns in rainfall–runoff. The hybrid model was useful in forecasting runoff.

In this study an attempt is made to model the variation of groundwater budget in the Azarshahr Plain aquifer using a conjugated wavelet-based artificial neural network (WANN)-GP model. This work differs from previously reported works in the sense that emphasis is given to improving the insight into the groundwater budget change and also linking the wavelet and GP for groundwater modelling.

## METHODOLOGY

### Study area

The Azarshahr Plain is one of the Urmia Lake sub-basins, and is located in Azerbaijan province, northwest Iran. The study area is densely populated, with 100% of its drinking, domestic and industrial water and 80% of agricultural water supplied from groundwater resources.

Azarshahr Chay is the main stream flow in the study area, which originates from Sahand Mountain and rarely discharges into the lake due to percolation and evaporation losses, as well as diversion of water for irrigation. The average annual precipitation of the study area is about 221.2 mm for the long-term period 1982 to 2009, whereas the annual evaporation is about 1,500 mm.

*et al.*2010). During the past decade, the study area has faced groundwater depletion and subsequently a reduction in reservoir volume, coinciding with population growth, great changes in climate and over-extraction of water resources in the study area, causing environmental hazards. Figure 2 shows the aquifer domain and piezometer locations in the study area.

### Wavelet transforms

As a contemporary tool of applied mathematics, wavelet transform (WT), is a signal processing strategy that has indicated higher performance contrasted with Fourier transform (FT) and short time FT in examining non-stationary signals. WT analysis, created during recent decades in the mathematics community, appears to be a more effective device than the FT in studying non-stationary time series (Partal & Kisi 2007). The principal point of preference for WTs is their capacity at the same time to get information on the time, location and frequency of a signal, while the FT will just give the frequency information of a signal.

The WT is implemented through discrete and continuous WT (DWT and CWT). Since different scales should be taken into consideration in CWT and using a numerical method, the equation integration for each scale is resolved. Calculating the wavelet coefficient is time-consuming in all scales and produces huge amounts of data. In other words, we can say that CWT consists of redundant and inefficient sections, which are its weak points (Adamowski 2007; Partal 2009); whereas the DWT has eradicated the CWT drawbacks. Meanwhile, it is an efficient alternative for the discrete data.

In DWT, the original time series is passed through high-pass and low-pass ﬁlters (digital filtering), getting time-scale signals. The results of digital filtering are detailed coefﬁcients and approximation series, obtained with the wavelet algorithm (Zhang & Li 2001). Every time that this procedure is repeated, the approximation and one or more details are gained.

*s*is the wavelet scale,

*t*is the time, τ is the translation parameter, while

*m*and

*n*are integers that control, respectively, the scale and time;

*s*

_{0}is a speciﬁed ﬁxed dilation step greater than 1; and τ

_{0}is the location parameter, which must be greater than zero. The mother wavelet is the transforming function. In the DWT, a time-scale signal is obtained using digital ﬁltering techniques.

Performing the above-mentioned transform, the raw data are divided into approximation (A) and details (D). The approximation consists of high scale and low frequency components of the signal. The details consist of low scale and high frequency components of the signal, which are obtained from low-pass and high-pass filters, respectively.

Consequently, the DWT was used to decompose the time series data belonging to the groundwater level for the wavelet analysis-artificial neural network (WA–ANN) models developed in this study.

### Hybrid WANN

The WT is appropriate in significant and potentially beneficial data mining, available in experimental sciences (prediction, reanalysis, global climate model simulations, etc.). Providing obvious information in a readable form, it can be applied to resolve analytical, classiﬁcation or forecasting issues. In a review of the applications of the WT in hydrologic time series modelling, Sang (2013) highlighted the complex information that can be drawn from such analysis: characterization and understanding of hydrologic series' multi-temporal scales, identiﬁcation of seasonality and trends, and data de-noising. Consequently, better interpreting of hydrological processes is derived from the decomposing ability of the WT (Nason & Sachs 1999; Adamowski 2008; Adamowski *et al.* 2009; Kisi 2010; Mirbagheri *et al.* 2010; Sang 2012).

Since AI has shown promise in modelling and forecasting non-linear hydrological processes and in handling large amounts of dynamicity and noise concealed in datasets, hybrid modelling of AI was employed for precise simulation of water level in piezometers and elimination of their statistical gaps in a long-term period. For this, the WT model was linked to ANN, producing WANN.

*i*is of the form: where

*w*is the weight vector of unit

_{ji}*i*and k is the number of neurons in the layer above the layer including unit

*i*.

*y*is the output from unit j, and

_{j}*y*is the bias of unit

_{i}*i*. This weighted sum Z, which is called the incoming signal of unit

*i*, is then passed through a transfer function

*f*to yield the estimates for unit

*i*.

In this study, DWT using Mallat's (1998) algorithm was used for decomposing the time series signal. The time series signal in this study is the water table fluctuation in piezometers, used only as the mother signal, which must be decomposed. The multi-resolution analysis by Mallat's algorithm generates approximations and details for a given time series signal. The general trend of the original signal and high frequency components are held and depicted by an approximation and detail, respectively. This results in breaking down the original signal into lower resolution constituents. Nourani *et al.* (2009) introduced the L = Int [log (N)] for choosing the number of decomposition levels or DWs, where L was the decomposition level while N was the number of time series data.

N-level DWT decomposes a signal x (t) into D_{1}, D_{2}…D_{N} and A_{N}, where D_{1} to D_{N} are details and A_{N} is an approximation. D_{1}, D_{2}…D_{N} and A_{N} are used as input to the ANN. The second step corresponds to training and testing phases using the ANN.

The WANN algorithm is summarized as follows:

Step 1: Multilevel wavelet analysis using DWT decomposes a signal into details (D

_{1}, D_{2}… D_{N}) and approximation (A_{N}), where N is the decomposition level. Water table time series data were decomposed into details D_{1}and D_{2}and an approximation A_{1}in this study. Decomposition levels have been selected with respect to the number of data used for each piezometer (the data of water table used for each piezometer) and they are shown in Table 1. For decomposition level, DL = log(No. Data) formula was used, following the suggestion of Wang & Ding (2003), Partal & Kisi (2007) and Nourani*et al.*(2009).Step 2: ANN is trained and tested using the details and approximation as input and the model performance is evaluated.

Piezometer | 3 | 4 | 5 | 6 | 7 | 8 | 10 | 11 | 12 | 14 | 19 | 21 | 23 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

No. data^{a} | 236 | 249 | 141 | 141 | 200 | 196 | 238 | 236 | 281 | 141 | 217 | 141 | 141 |

Wavelet type | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 |

Level | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |

Hidden layer | 8 | 8 | 6 | 6 | 8 | 8 | 8 | 8 | 8 | 6 | 8 | 6 | 6 |

Piezometer | 3 | 4 | 5 | 6 | 7 | 8 | 10 | 11 | 12 | 14 | 19 | 21 | 23 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

No. data^{a} | 236 | 249 | 141 | 141 | 200 | 196 | 238 | 236 | 281 | 141 | 217 | 141 | 141 |

Wavelet type | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 | db4 |

Level | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 |

Hidden layer | 8 | 8 | 6 | 6 | 8 | 8 | 8 | 8 | 8 | 6 | 8 | 6 | 6 |

^{a}The number of data for each piezometer (number of months).

### Genetic programming

GP is a kind of artiﬁcial intelligence method that is based on the random iterative searching process to achieve an appropriate relationship between input and output. The common structure of this method is the tree shape, representing the expression. Variables, functions and operators in this structure are situated in the nodes, which are linked together by branches.

GP is an evolutionary algorithm based on Darwinian theories of natural selection and survival of the fittest. The algorithm considers an initial population of randomly generated equations, derived from the random combination of input variables, random numbers and functions. The function can include arithmetic operators (plus, minus, multiply and divide), mathematical functions (sin, cos, exp, log) etc., which have to be chosen based on some understanding of the process. This population is then subjected to an evolutionary process and the fitness of the evolved programs are evaluated; individual programs that best fit the data are then selected from the initial population. The programs that best fit are selected to exchange part of the information between them to produce a better program through ‘crossover’ and ‘mutation’. The user must decide a number of GP parameters before applying the algorithm to the data. The program that fitted the data less well is discarded. This evolution process is repeated over successive generations and is driven towards finding symbolic expressions describing the data, which can be scientifically interpreted to derive knowledge about the process being modelled (Sivapragasam *et al.* 2015).

However, in this study, the simulated water tables by WANN were used as inputs to the GP model for time series prediction as can be seen from the flowchart of the study steps (Figure 1). It means that the simulated water tables of each piezometer have been entered into the GP model and forecasting has been done after that. For this aim, GeneXproTools was utilized and the training, testing and forecasting of 12-month ahead water tables was done.

### Performance evaluation

^{2}), RMSE, MAE and Nash–Sutcliffe efficiency (NSE). The mathematical form of performance criteria are as below: where and are observation and forecasted values, respectively, and n is the number of data.

### Groundwater budget calculation

The procedure of monthly water level forecasting is as below:

Forecasted water level for each month = [[(polygon area of 1st piezometer × 1st piezometer water level) + (polygon area of 2nd piezometer × 2nd piezometer water level) + ··· + (polygon area of Nth piezometer × Nth piezometer water level)]/(polygon area summation)]

In hydrological studies, we need to know the water budget of the groundwater reservoir, where sometimes we do not have the data of the input and output parameters of the study area, such as precipitation, evaporation, etc.; in this situation, we can evaluate the reservoir groundwater budget by calculating the fluctuations in the water table of the aquifer during the given period (here from 2012 to 2013). Knowing the difference between the water table at the beginning and end of the period, and also knowing the specific storage, S, and the area of the reservoir (A), it is possible to calculate the changes of reservoir groundwater volume. This leads us to understand the groundwater budget, i.e., if the reservoir gained or lost water during the given period of time.

*h*is the water table change during the study period, forecasted by the hybrid GP-WANN model in this study.

## RESULT AND DISCUSSION

*et al.*(2009) and considering the statistical period (10 to 20 years) of the water level in the piezometers, two decomposition levels were used in DWT. Thus, the original groundwater budget time series were decomposed into D1, D2 and A2, where D1 and D2 are details and A2, is an approximation. D1, D2 and A2 were inputted into the ANN. The number of data used for ANN varies from 141 to 249, including water levels for a 10–20 year period, and 80% and 20% of data were selected for training and testing the ANN, respectively. The Levenberg–Marquardt algorithm was chosen as the training algorithm. Different hidden node numbers were tried, and the optimal hidden nodes were found to vary between 2 and 8 for the optimal ANN models. Simulation iterations were finished at 23 to 25 to achieve the best results. Figure 6 depicts the testing results for 13 piezometers in the aquifer domain. It is noteworthy that piezometer numbers 4, 8 and 14 had 1, 1 and 3 month gaps during the test data period. WANN covered the gaps and forecasted the water levels in the test period.

The evaluation of model performance was done by statistical criteria (R^{2}, RMSE, MAE and NSE) and the results are shown in Table 2.

Piezometer | 3 | 4 | 5 | 6 | 7 | 8 | 10 | 11 | 12 | 14 | 19 | 21 | 23 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

R^{2} | 0.97 | 0.92 | 0.9 | 0.88 | 0.9 | 0.81 | 0.87 | 0.83 | 0.81 | 0.96 | 0.94 | 0.8 | 0.84 |

RMSE | 0.0698 | 0.2365 | 0.1036 | 0.7104 | 0.0759 | 4.45 | 0.3869 | 0.556 | 0.622 | 0.2683 | 0.3578 | 1.1643 | 0.4605 |

MAE | 0.0493 | 0.1603 | 0.0691 | 0.0156 | 0.538 | 3.0362 | 0.2896 | 0.3633 | 0.4395 | 0.2095 | 0.1659 | 0.686 | 0.2922 |

NSE | 0.93 | 0.91 | 0.92 | 0.83 | 0.91 | 0.75 | 0.78 | 0.85 | 0.82 | 0.96 | 0.94 | 0.81 | 0.74 |

Piezometer | 3 | 4 | 5 | 6 | 7 | 8 | 10 | 11 | 12 | 14 | 19 | 21 | 23 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

R^{2} | 0.97 | 0.92 | 0.9 | 0.88 | 0.9 | 0.81 | 0.87 | 0.83 | 0.81 | 0.96 | 0.94 | 0.8 | 0.84 |

RMSE | 0.0698 | 0.2365 | 0.1036 | 0.7104 | 0.0759 | 4.45 | 0.3869 | 0.556 | 0.622 | 0.2683 | 0.3578 | 1.1643 | 0.4605 |

MAE | 0.0493 | 0.1603 | 0.0691 | 0.0156 | 0.538 | 3.0362 | 0.2896 | 0.3633 | 0.4395 | 0.2095 | 0.1659 | 0.686 | 0.2922 |

NSE | 0.93 | 0.91 | 0.92 | 0.83 | 0.91 | 0.75 | 0.78 | 0.85 | 0.82 | 0.96 | 0.94 | 0.81 | 0.74 |

It can be clearly seen that the greater the fluctuation of the water level, the lower the performance of the model. The lower performance belongs to piezometer numbers 8, 11, 12 and 21, and may be driven by their situation. For example, piezometer No. 8 is located near the road and belongs to a tubing manufacturer that uses groundwater for its production, and therefore the water level experiences oscillation. Also piezometers numbers 11 and 12 are near the river and perhaps affected by river fluctuations. Thus it can be inferred that the model performs with lower accuracy and more error. Table 3 also shows the results of the predicted water level for the 13 piezometers and their Thiessen polygon area used for groundwater budget computing.

No.^{a} | Polygon area | H1^{b} | H2 | H3 | H4 | H5 | H6 | H7 | H8 | H9 | H10 | H11 | H12 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

3 | 437,212.86 | 1,332.04 | 1,332.1 | 1,332.15 | 1,332.16 | 1,332.12 | 1,332.11 | 1,332.08 | 1,331.97 | 1,331.96 | 1,331.93 | 1,331.94 | 1,331.93 |

4 | 2,143,656.4 | 1,291.8 | 1,291.75 | 1,291.95 | 1,292.18 | 1,292.35 | 1,292.42 | 1,292.33 | 1,292.45 | 1,292.58 | 1,292.28 | 1,292.25 | 1,292.18 |

5 | 4,218,328.7 | 1,321.21 | 1,321.21 | 1,321.4 | 1,321.51 | 1,321.53 | 1,321.58 | 1,321.59 | 1,321.69 | 1,321.7 | 1,321.63 | 1,321.62 | 1,321.61 |

6 | 2,943,760.8 | 1,281.19 | 1,281.06 | 1,281.17 | 1,281.15 | 1,281.02 | 1,280.98 | 1,280.99 | 1,281.11 | 1,281.13 | 1,281.09 | 1,281.12 | 1,281.05 |

7 | 11,785,231 | 1,272.09 | 1,271.68 | 1,272.58 | 1,274.7 | 1,275.51 | 1,275.21 | 1,275.18 | 1,273.75 | 1,272.97 | 1,272.81 | 1,271.8 | 1,271.7 |

8 | 3,642,347.3 | 1,357.56 | 1,354.21 | 1,349.95 | 1,349.48 | 1,351.11 | 1,351.28 | 1,357.96 | 1,358.41 | 1,358.39 | 1,358.37 | 1,358.42 | 1,357.86 |

10 | 7,371,027.7 | 1,437.31 | 1,437.35 | 1,438.61 | 1,438.87 | 1,439.00 | 1,438.12 | 1,438.45 | 1,438.78 | 1,438.73 | 1,438.61 | 1,437.05 | 1,436.63 |

11 | 11,164,588 | 1,334.75 | 1,334.65 | 1,334.9 | 1,335.4 | 1,336.12 | 1,336.92 | 1,337.36 | 1,337.44 | 1,337.22 | 1,337 | 1,335.51 | 1,334.92 |

12 | 9,018,546.3 | 1,385.45 | 1,386.59 | 1,386.27 | 1,387.96 | 1,388.85 | 1,389.85 | 1,390.59 | 1,391.42 | 1,391.32 | 1,391.19 | 1,386.03 | 1,385.79 |

14 | 7,413,532 | 1,283.02 | 1,283.24 | 1,283.89 | 1,284.57 | 1,285.14 | 1,285.62 | 1,285.27 | 1,285.11 | 1,284.84 | 1,284.76 | 1,284.17 | 1,283.94 |

19 | 6,827,615.3 | 1,278.25 | 1,278.23 | 1,278.22 | 1,278.2 | 1,277.57 | 1,275.27 | 1,275.35 | 1,275.39 | 1,275.38 | 1,275.33 | 1,275.27 | 1,275.21 |

21 | 5,162,200.5 | 1,416.43 | 1,415.73 | 1,415.45 | 1,415.68 | 1,415.83 | 1,416.89 | 1,416.91 | 1,416.7 | 1,416.63 | 1,416.57 | 1,415.84 | 1,416.53 |

23 | 9,524,845.5 | 1,275.44 | 1,275.32 | 1,275.26 | 1,275.97 | 1,276.26 | 1,276.77 | 1,277.23 | 1,277.00 | 1,276.9 | 1,276.86 | 1,276.09 | 1,275.97 |

No.^{a} | Polygon area | H1^{b} | H2 | H3 | H4 | H5 | H6 | H7 | H8 | H9 | H10 | H11 | H12 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

3 | 437,212.86 | 1,332.04 | 1,332.1 | 1,332.15 | 1,332.16 | 1,332.12 | 1,332.11 | 1,332.08 | 1,331.97 | 1,331.96 | 1,331.93 | 1,331.94 | 1,331.93 |

4 | 2,143,656.4 | 1,291.8 | 1,291.75 | 1,291.95 | 1,292.18 | 1,292.35 | 1,292.42 | 1,292.33 | 1,292.45 | 1,292.58 | 1,292.28 | 1,292.25 | 1,292.18 |

5 | 4,218,328.7 | 1,321.21 | 1,321.21 | 1,321.4 | 1,321.51 | 1,321.53 | 1,321.58 | 1,321.59 | 1,321.69 | 1,321.7 | 1,321.63 | 1,321.62 | 1,321.61 |

6 | 2,943,760.8 | 1,281.19 | 1,281.06 | 1,281.17 | 1,281.15 | 1,281.02 | 1,280.98 | 1,280.99 | 1,281.11 | 1,281.13 | 1,281.09 | 1,281.12 | 1,281.05 |

7 | 11,785,231 | 1,272.09 | 1,271.68 | 1,272.58 | 1,274.7 | 1,275.51 | 1,275.21 | 1,275.18 | 1,273.75 | 1,272.97 | 1,272.81 | 1,271.8 | 1,271.7 |

8 | 3,642,347.3 | 1,357.56 | 1,354.21 | 1,349.95 | 1,349.48 | 1,351.11 | 1,351.28 | 1,357.96 | 1,358.41 | 1,358.39 | 1,358.37 | 1,358.42 | 1,357.86 |

10 | 7,371,027.7 | 1,437.31 | 1,437.35 | 1,438.61 | 1,438.87 | 1,439.00 | 1,438.12 | 1,438.45 | 1,438.78 | 1,438.73 | 1,438.61 | 1,437.05 | 1,436.63 |

11 | 11,164,588 | 1,334.75 | 1,334.65 | 1,334.9 | 1,335.4 | 1,336.12 | 1,336.92 | 1,337.36 | 1,337.44 | 1,337.22 | 1,337 | 1,335.51 | 1,334.92 |

12 | 9,018,546.3 | 1,385.45 | 1,386.59 | 1,386.27 | 1,387.96 | 1,388.85 | 1,389.85 | 1,390.59 | 1,391.42 | 1,391.32 | 1,391.19 | 1,386.03 | 1,385.79 |

14 | 7,413,532 | 1,283.02 | 1,283.24 | 1,283.89 | 1,284.57 | 1,285.14 | 1,285.62 | 1,285.27 | 1,285.11 | 1,284.84 | 1,284.76 | 1,284.17 | 1,283.94 |

19 | 6,827,615.3 | 1,278.25 | 1,278.23 | 1,278.22 | 1,278.2 | 1,277.57 | 1,275.27 | 1,275.35 | 1,275.39 | 1,275.38 | 1,275.33 | 1,275.27 | 1,275.21 |

21 | 5,162,200.5 | 1,416.43 | 1,415.73 | 1,415.45 | 1,415.68 | 1,415.83 | 1,416.89 | 1,416.91 | 1,416.7 | 1,416.63 | 1,416.57 | 1,415.84 | 1,416.53 |

23 | 9,524,845.5 | 1,275.44 | 1,275.32 | 1,275.26 | 1,275.97 | 1,276.26 | 1,276.77 | 1,277.23 | 1,277.00 | 1,276.9 | 1,276.86 | 1,276.09 | 1,275.97 |

^{a}No. = piezometer number.

^{b}H1…H12 = water level from 1st month to 12th month.

Month | Sep-12 | Oct-12 | Nov-12 | Dec-12 | Jan-13 | Feb-13 | Mar-13 | Apr-13 | May-13 | Jun-13 | Jul-13 | Aug-13 |
---|---|---|---|---|---|---|---|---|---|---|---|---|

Water level | 1,326.63 | 1,326.49 | 1,326.60 | 1,327.33 | 1,327.77 | 1,327.85 | 1,328.35 | 1,328.25 | 1,328.06 | 1,327.95 | 1,326.69 | 1,326.51 |

Month | Sep-12 | Oct-12 | Nov-12 | Dec-12 | Jan-13 | Feb-13 | Mar-13 | Apr-13 | May-13 | Jun-13 | Jul-13 | Aug-13 |
---|---|---|---|---|---|---|---|---|---|---|---|---|

Water level | 1,326.63 | 1,326.49 | 1,326.60 | 1,327.33 | 1,327.77 | 1,327.85 | 1,328.35 | 1,328.25 | 1,328.06 | 1,327.95 | 1,326.69 | 1,326.51 |

Using Equation (7), the predicted water table for the water budget domain, considering the aquifer domain area (81.65 km^{2}) and storage coefficient, derived about 0.036 from pumping wells and qanat discharges in the Azarshahr Plain detailed collected data, done by the Azerbaijan Territorial Water Association (ATWA) in 2009. The reservoir volume, the groundwater budget of the study area, has reduced by about 0.3557 million cubic metres (MCM) during the prediction period (September 2012 to August 2013).

## CONCLUSION

This paper has given an account for the widespread use of AI in the simulation and prediction of environmental processes. Hybrid modelling of wavelet and ANN in simulating the water table and then forecasting the future water table using GP was applied to determine the study area groundwater reservoir changes during the forecasted period of time. Performance assessment shows satisfying results, revealing the accuracy of AI hybrid models. Not only the WANN hybrid model but also the GP model has shown its ability in simulation and prediction of natural events. The next year's water budget was measured, using the Thiessen polygon method, knowing the area of the aquifer domain and its storage coefficient and forecasting the aquifer water table fluctuation. The groundwater reservoir has lost about 0.35 MCM of its storage during the predicted period, which indicates the importance of monitoring the groundwater resource and, of course, the predicted water budget can be taken into account for future environmental plans.