Evaluating the Impacts of Pumping on Aquifer Depletion in Arid Regions Using MODFLOW, ANFIS and ANN

Abstract

In arid regions, the groundwater drawdown consistently increases, and even for a constant pumping rate, long-term predictions remain a challenge. The present research applies the modular three-dimensional finite-difference groundwater flow (MODFLOW) model to a unique aquifer facing challenges of undefined boundary conditions. Artificial neural networks (ANN) and adaptive neuro fuzzy inference systems (ANFIS) have also been investigated for predicting groundwater levels in the aquifer. A framework is developed for evaluating the impact of various scenarios of groundwater pumping on aquifer depletion. A new code in MATLAB was written for predictions of aquifer depletion using ANN/ANFIS. The geotechnical, meteorological, and hydrological data, including discharge and groundwater levels from 1980 to 2018 for wells in Qassim, were collected from the ministry concerned. The Nash–Sutcliffe efficiency and mean square error examined the performance of the models. The study found that the existing pumping rates can result in an alarming drawdown of 105 m in the next 50 years. Appropriate water conservation strategies for maintaining the existing pumping rate can reduce the impact on aquifer depletion by 33%.

1. Introduction

Groundwater modeling for sustainable management is an important field in regions that use groundwater as their main water supply. The matter becomes more serious where the recharge to the aquifer is limited and the groundwater is depleting at an alarming rate (Figure 1). So, it is a hot issue and relevant research has been addressing various issues of groundwater worldwide [1,2,3]. However, there are several challenges in groundwater modeling used to provide some of the basic information for engineering designs, sustainable groundwater planning, and environmental impact assessment. Challenges exist in understanding system dynamics and the flow patterns, analyzing the groundwater responses to stresses in the system, and assessing the sustainable groundwater yield. In short, groundwater modeling is an area of continued research, as identified by several publications [4,5,6,7,8,9]. The groundwater flow is governed by the well-known Boussinesq equation. A variety of computer programs have been developed to solve this equation. Mainly, finite element or finite difference techniques are used for this purpose. MODFLOW, PLASM, and AQUIFEM-1 to AQUIFEM-N are a few examples of such software [8,10]. MODFLOW provides the most commonly used solution of the equation based on finite differences. This model has been used by several researchers from the early 1980s to 2020 [4,11,12,13,14,15]. In spite of the long-term use of numerical modeling for groundwater, there remains space for further research in this field of specialization because of the unique problems with various aquifers [16,17].

The modeling of aquifers becomes highly challenging if the boundary conditions are not well defined for their certain parts. Pujades et al. [18,19] comprehensively identified the difficulties associated with solving the problems of groundwater flow modeling with unknown boundary conditions and aquifer parameters. They have sought solutions to these problems by using approximate and numerical methods applicable to a specific groundwater flow situation characterized by its hydraulic behavior of diaphragm walls of cut and cover tunneling and an underground enclosure. They used pumping tests data to obtain the aquifer parameters (hydraulic conductance and specific storage) by the Theis method [20]. Pujades et al. [21] used pumping tests data to obtain the same aquifer parameters by Cooper and Jacob’s method [22]. Although the previous studies have provided guidelines to solve the problem of unknown boundary conditions and for the estimation of aquifer parameters, their pragmatism can only be evaluated by applying them to diverse cases with site-specific constraints [23,24,25,26].

A part of the Saq Aquifer in Qassim, Saudi Arabia also has no well-defined boundary conditions. In the first objective of the present study, we addressed the problem of the boundary conditions and other aquifer parameters using hit and trial and the calibration subroutine of MODFLOW. There is also a need to address the above-mentioned problem through alternative, relatively simpler, modeling techniques. Several artificial intelligences based techniques have been found as useful alternatives, e.g., artificial neural networks (ANN), adaptive neuro fuzzy inference systems (ANFIS), particle swarm optimization-neuro fuzzy inference systems (PSO-ANFIS), support vector regression (SVR), bidirectional long short-term memory (Bi-LSTM) network, and gated recurrent unit (GRU). Elbaz et al. [27,28] integrated PSO-ANFIS with genetic algorithm models to predict the earth pressure balance for the performance of a shield via tunneling at a site in China. The very small error found between the observed and simulated earth pressure balance showed the high performance of their models. Recently, Gao et al. [29] used GRU for solving a similar problem in Luoyang, China. Although the application of these deep learning arterial intelligence models in the field of groundwater predictions is relatively rare to date, some researchers used the ANN, ANFIS, and SVR to simulate groundwater-levels [30,31,32,33,34,35,36].

Panahi et al. [37] used convolutional neural network and SVR models to map groundwater but did not address the challenges associated with the in-depth-analysis of groundwater levels. Fan et al. [38] applied machine learning methods to analyze pumping tests to determine the aquifer parameters. They concluded that the machine learning models can successfully help to analyze the pumping tests for the estimation of aquifer parameters. Chen et al. [39] compared the application of machine learning in groundwater modeling with numerical techniques. They found that the performance of machine learning models (ANN and SVR) is better than that of hydraulic models for short-term predictions. However, hydraulic techniques of groundwater simulations are based on physical laws and these are much better methods than the artificial intelligence models in general (

Table 1. Comparison of GCM with data-driven models ANN and ANFIS.

Characteristics	Hydraulic Model	ANN/ANFIS	Remarks
Computational expenses (need high speed computers and time for computations)	High	Moderate	- Comparatively higher computing facilities are required for hydraulic models. - Moderate computing power may function adequately for the ANN. [8,39]
Data and model processing requirements	Complex	Moderately complex	- Requires precise values of aquifer parameters and high resolution mesh with lower computational time steps. - Data of a specific aquifer are used to predict the future water levels and quality parameters in ANN/ANFIS application. [8,10,32,33,39]
Model Type	Physically distributed	Black-box (data driven)	- Hydraulic models incorporate physical processes and laws involved in predicting groundwater levels and quality parameters. - ANN/ANFIS are data-driven models and do not involve any physical processes. These models only use the data observed during training, testing, and predictions. [10,32,33,39]
Model bias	Moderate to good	Good for specific data but limited in generality	- The biases are of varying magnitude (wide range of variation) in hydraulic model projections. - The biases are limited to a certain range in most cases of data-driven model projections, but ANNs/ANFIS lack generality. [8,10,32,33,39]
Long-term future predictions	Easy, once the model is calibrated	Challenging	- Long-term future predictions are easy once the hydraulic model is calibrated and validated. - High expertise and experience of specific problem is required for long-term future predictions by ANN/ANFIS. Expertise is required for the normalization and de-normalization of data. [8,10,32,33,39]

). Nevertheless, extensive aquifer data is required for use of these models. Their study showed that further research is required to further explore the use of ANN in the field of groundwater modeling. Accordingly, the second objective of the present paper is to find a simple solution to the problem of undefined boundary conditions of aquifers using ANN and ANFIS models.

The third objective of this study is the use of ANN/ANFIS models for long-term future predictions of groundwater levels to explore the aquifer depletion. Most of the Gulf countries are water-stressed and require sustainable groundwater planning and management [40]. The Saq Aquifer in the Qassim Region of Saudi Arabia represents a unique characteristic of arid regions, i.e., the aquifer recharge is almost negligible, but the groundwater pumping is continuous to meet municipal and irrigation requirements. The entire aquifer can ultimately dry sometime in the future if the present pumping trends continue without appropriate remedial measures. A simple but reliable groundwater model is required for this purpose. Hence, three ANN and an ANFIS models were investigated in this study and the best one was used for long-term predictions of groundwater-levels. It is worth mentioning here that most of the modelling based on ANN/ANFIS in past research dealt with only the short-term predictions of groundwater-levels [34,35,36,37,41,42]. Khedri et al. [42]. applied five various data-driven techniques for the simulation of groundwater levels. Their study only dealt with the short-term predictions of groundwater levels. According to them, models based on artificial intelligence could predict groundwater levels properly for one to two future months, but in the case of forecasting for three future months, the performance of these techniques was substandard. Long-term predictions using ANN/ANFIS is a challenging and novel research area, which is addressed properly in this paper. A state-of-the-art model was developed by writing a code in MATLAB for the long-term simulation of groundwater-levels using ANN/ANFIS.

Although there is no conceptual difference between our model and the classical ANN model, in the classical neural network toolbox, the user cannot, (i) define initial weights of the links and bias values, (ii) define different types of activation functions between the layers, (iii) outline different training algorithms (the 12 algorithms), and (iv), most importantly, users don’t have the option to select the type of the division of the database for the cross validation into subsets. Our new code facilitates defining the initial weights of the links and bias values, changing the number of hidden layers, neurons in each hidden layer, activation functions, training algorithms, and division of data base (i.e., randomly, index, group) for the cross validation into subsets (training, validation and testing). The new ANN code is developed in a MATLAB environment due to the limitations of neural network toolbox. As author used different combinations of the above-mentioned parameters, the results of the present study show an optimized architecture of the ANN model.

It is worth mentioning here that long-term predictions for time series problems no doubt require deep learning models, such as recurrent neural networks (RNN), LSTM, convolutional neural networks (CNN), etc. For long-term predictions, we adopted a simple approach by changing the targets and predictions into groundwater-level-changes, instead of groundwater-levels, and make it a conventional ANN/ANFIS simulation problem. Hence, we have tackled the problem of long-term simulations by training the ANN/ANFIS model for the prediction of changes in groundwater levels instead of the direct simulation of water levels. The normalization and de-normalization were made for target and predicted groundwater level changes, respectively. This makes the problem of long-term predictions similar to that dealt with the classical ANN/ANFIS models for short-term predictions. The impacts of various pumping rates on groundwater levels are studied to find an optimal solution for reducing the depletion of Saq aquifer.

The results of ANFIS are added in this research to check the efficiency of a hybrid type artificial intelligence model. ANFIS integrates the neural networks with the fuzzy rule-based system to accommodate the system uncertainties. First, this integrated approach develops a fuzzy model that derives its input variables from the fuzzy rules defined based on the input-output data. Next, the neural network tweaks these rules and generates the final ANFIS model.

2. Materials and Methods

2.1. Methodological Framework

In the present study, we used two types of techniques, including the hydraulic model MODFLOW and the data-driven models (ANN and ANFIS), to explore the Saq Aquifer depletion. A brief comparison of these models is given in Table 1. A framework illustrating the overall methodology is shown in Figure 2. The method used in this study comprised three main steps. The first step consisted of data collection regarding groundwater-levels, pumping rates, well locations, and aquifer parameters of the study area. We normalized the data regarding pumping rates and groundwater-levels for ANN and ANFIS models. The second step was to calibrate/train/test/validate the hydraulic and ANN/ANFIS models. Part of the data (say 23 years out of 38 years, i.e., 60%) was used for training/calibration and the remainder was used for testing/validation. So, the data were divided into three parts, namely 60%, 20%, and 20%, respectively, for the training, testing, and validation of ANN/ANFIS models and two parts, 60% and 40%, for the calibration and validation of the hydraulic model. Coding was prepared in MATLAB to develop an ANN and ANFIS to accomplish the required objectives. In the third step, the models were applied for future predictions of groundwater levels to study the aquifer depletion. The rate of groundwater use is increasing day after day because of the population increase and the ever-rising standards of living. According to a report from the Ministry of Water and Electricity Saudi Arabia, the population growth rate is about 3% in the Riyadh and Qassim Regions [43]. A total of four future scenarios for various pumping rates were examined to explore the aquifer depletion. Three scenarios were built with an increase in the prevailing water extraction rate by 1%, 2% and 3% per year from 2020 to 2070. The fourth scenario was examined by continuing a constant rate of pumping from 2018 to 2070 without any change in the pumping rate as observed in 2018.

2.2. Study Area

Our study area included one of the most important cities of the Al-Qassim region, the Buraydah, and its surroundings, as shown in Figure 3. Buraydah is the capital of the region and is located between longitudes 43°50′55.99″ and 44° 8′ 59.85″ E and latitudes 26°10′47.28″ N and 26°27′5.60″ N. It has mainly dry weather [44]. The average annual rainfall in the area is about 125 mm [4,45].

The water supply to the region is mainly dependent upon groundwater. The Saq Aquifer, which is one of the well-known aquifers in Saudi Arabia, is the major source of water supply. The cross-section of the Saq Aquifer is shown in Figure 4. The outcrop of the Saq aquifer is very vast. It encompasses nearly 1200 km only in Saudi Arabia, joining the Jordanian border as far south as longitude 45° E and latitude 24°30′ N. It is a confined aquifer with a variable thickness changing from 700 m in its northern part to 400 m in its southern part. The average thickness of the aquifer in the study area is nearly 500 m [4,45].

The subsurface strata comprise medium- to course-sized sandstone. There are some local areas containing fine sandstone also. The Saq aquifer groundwater has an electrical conductivity of approximately 0.276 S/m. The pumping wells are partially penetrating, having an average depth of about 650 m. The screen of the wells is 125 m long [4].

2.3. Description of Hydraulic Model (MODFLOW)

The well-known Bossinesq equation governs groundwater flow is expressed as:

$$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)=S_s\frac{\partial h}{\partial t}-W$$

(1)

where the values of the hydraulic conductivities are represented by Kxx, Kyy, and Kzz in the x, y, and z directions, respectively, h represents the groundwater-head, the specific storage of the aquifer is expressed by Ss, W represents a source/sink, and time is expressed by t. The finite difference scheme is applied to obtain the numerical solution to the above equation in the MODFLOW. The input and output data for MODFLOW is presented in

Table 2. Input and output data for MODFLOW.

Data Name	Description	Value
Input Data
Hydraulic conductivity	Estimated by pumping test and Cooper Jacob Method [4,22]	Kxx = Kyy = Kzz = 3.6 m/day
Transmissivity	Estimated by pumping test and Cooper Jacob Method [4,22]	4500 m2/day
No of aquifer layers	Estimated by previous studies [4]	5 (631, 125, 125, 125, and 125 m)
Mesh size	Adopted from previous studies [4]	30 km × 20 km area divided into 80 rows × 70 columns
Pumping rates	From Ministry of Environment, Agriculture and Water, Saudi Arabia	For 55 pumping/observational wells (1980 to 2018)
Water levels (a)	From Ministry of Environment, Agriculture and Water, Saudi Arabia	For 55 pumping/observational wells (1980 to 2018)
Initial conditions	Adopted from previous studies [4]	Water levels at the start time of simulations
Boundary conditions	Calibrated assuming general head boundary conditions	Challenging task
Output
Water levels (b)	Predicted by the model	For 55 pumping/observational wells (1980 to 2018) for calibration and validation
Water levels (c)	Predicted by the model	For 55 pumping/observational wells (1918 to 2070) future predictions

. Selection of the mesh size and appropriate boundary conditions is highly challenging task in the case of the Saq aquifer. In this paper, a rectangular type mesh of 30 × 20 km was developed by creating 80 rows and 70 columns. The number of rows and columns in the mesh were increased significantly to obtain denser mesh near the pumping wells for simulating the groundwater levels as precisely as possible. As mentioned earlier, the aquifer thickness (in the vertical direction) was 500 m, above which there was an impermeable cover. The upper layer above the 500 m thick aquifer is about 628 m [4,45]. The average ground surface level is around 600 m from the mean sea level. The pumping wells are partially penetrating wells with a screen length of 125 m, so the confined part of the 500-m thick aquifer is divided into four layers, each of 125 m, to simulate the screen effects precisely. In this way, five layers with thicknesses of 631, 125, 125, 125, and 125 m were considered, while the top layer was considered impermeable (Figure 5).

The parameters of the aquifer below the 600 m soil cover were found by the pumping test. The transmissivity of the aquifer was found to be 4500 m²/day. The aquifer is confined, homogeneous, and isotropic. The hydraulic conductivity of the aquifer was found to be 3.6 m/day and was the same in all the three directions (Kxx = Kyy = Kzz). The recharge from the surface to the aquifer is negligible. The seepage from the basement and the deep groundwater outflow were considered negligible. The aquifer is very large and has known boundaries only at its ends and outcrops. There is continued pumping with negligible recharge, so the boundary conditions change continually. Hence, the choice of boundary conditions for the study area was a difficult task. The general head boundary condition was the appropriate choice and was adjusted by calibration. This is further explained in the calibration section.

2.4. Artificial Neural Networks (ANN)

The artificial neural network (ANN) model is a data-driven model [41,46]. The ANN is one of the artificial intelligence techniques that simulates the actions of the human brain with the help of neurons. ANNs are becoming very common for hydrologic modeling and have been used to solve problems of applied nature in engineering and science [47,48]. These models may be congregated into two key categories: feed-forward and feed-back networks [46,47,49]. The most commonly used family of feed-forward networks contains a layered-network in which the neurons are organized in layers with unidirectional connections between the layers. This is known as multilayer perceptron (MLP) [47,48,49]. The input and output data for ANN are presented in

Table 3. Input and output data/parameters for ANN.

Parameter Name	Description	Value/Comment
Input Data
Number of hidden layers	Adopted to investigate the best architecture of the model	Two to three
Number hidden neurons in each layer	Adopted to investigate the best architecture of the model	5 to 15 (5, 10, 15)
Activation functions	Tested various functions	-
Training algorithms	Tested various training algorithms	-
Data base division for training, testing and cross validation	Tested various divisions	-
Pumping rates	From Ministry of Environment, Agriculture and Water, Saudi Arabia	For selected pumping/observational wells (1980 to 2018)
Water levels (a)	From Ministry of Environment, Agriculture and Water, Saudi Arabia	For selected pumping/observational wells (1980 to 2018)
Output
Water levels (b)	Predicted by the model	For selected pumping/observational wells (1980 to 2018) for training/testing/validation
Water levels (c)	Predicted by the model	For 55 pumping/observational wells (1918 to 2070) future predictions

. Figure 6a illustrates an MLP with three types of layers, including input, output, and hidden layer. The neurons in the input layer act as buffers for allocating the input signals to neurons in the hidden layer. The input signals can be expressed as x_i (i = 1, 2, …, n with n being the total number of input signals). Each neuron j in the hidden layer characterizes the sum of its input signals x_i after weighting them with the strengths of the respective connections w_ij from the input layer. It then calculates its output y_j as a function f of the sum of w_ij × x_i.

Three of the built-in functions available in MATLAB (Levenberg Marquardt (LM), Bayesian regularization, and scaled conjugate gradient) were tested. The LM algorithm combines the merits of two training algorithms, namely the steepest descent and Gaussian–Newton methods, and searches for the global minima function to optimize the solution [41,47,48,49]. Bayesian regularization is a highly robust training algorithm for simulating the short-term fluctuations in groundwater levels. However, the Levenberg–Marquardt algorithm yields convergence comparatively more quickly than Bayesian regularization [49]. The scaled conjugate gradient is an efficient training algorithm with quick convergence and a high degree of accuracy [49,50,51].

The performance of the developed models was tested using Nash–Sutcliffe model efficiency (NSE) and mean square error (MSE), as explained in the coming sections. Finally, the results were generated for future groundwater levels using an ANN for Buraydah Qassim up to the year 2070 for various pumping rate scenarios.

2.5. Adaptive Neuro Fuzzy Inference System

ANFIS is a very popular branch of artificial intelligence models. A combination of fuzzy-if-then rules are used to describe the input-output relation of a real system in ANFIS [32,33]. The output of each of these rules is defined by a set of input variables and a constant term. The average of the weighted output of each rule becomes the concluding output. The functional block of ANFIS, consisting of fuzzification, inference, and defuzzification processes, is presented in Figure 6b [30,31].

2.6. Model Performance Evaluation

There are numerous statistical indices that are used for the evaluation of model performance. These indices are based on the comparison of measured and predicted data. It is worth noting that about 38 years (1980–2018) of data recorded for groundwater levels and pumping rates were used to check the performance of the models for various stages, including training/calibration, testing, and validation. A similar approach was adopted in previous studies [44,52].

The parameters adopted in this study to measure the performance of the models are given by the following equations [44,52,53]. The Nash–Sutcliffe model efficiency (NSE) is given as:

$$NSE=1-\frac{{\displaystyle {\sum }_{i=1}^n(}GW{L}_i^o-GW{L}_i^p{)}^2}{{\displaystyle {\sum }_{i=1}^n(}GW{L}_i^o-GW{L}_{avg}{)}^2}$$

(2)

where GWL is the variable representing the groundwater level, o is the observed value of GWL, p represents the predicted value, and n denotes the total number of data points. As per the criteria adopted by Rauf and Ghumman [54], values from 0.75 to 1.0 of NSE can be categorized as ‘very good’, 0.65 to 0.75 can be considered ‘good’, 0.5 to 0.65 as ‘satisfactory’, and the values between 0.4 and 0.5 represent an ‘acceptable’ performance of the model.

The mean bias error (MBE) is given as:

$$MBE={{\displaystyle \sum }}_{i=1}^n\frac{\left(GW{L}_i^p-GW{L}_i^o\right)}{n}$$

(3)

The positive values of MBE represent overestimated predictions whereas negative values are for underestimated predictions [54].

3. Results and Discussion

3.1. Hydraulic Model Results

The most important and challenging parameter in the calibration of the hydraulic model was the adjustment of the general head boundary condition. The results of calibration of MODFLOW are shown in Figure 7a,b. The NSE value was found to be 0.9. This value of NSE lies in the range of “very good” according to the criteria adopted by Rauf and Ghumman [54]. Furthermore, our results are similar to those of Mohanty et al. [55] and Lyons et al. [16], although the hydraulic conductivity and specific storage were the most sensitive parameters in the calibration phase of their groundwater hydraulic model, but adjustment of the boundary condition was most crucial issue in our case. The resultant predicted vs. observed line was close to the 45-degree line, which also showed that the results of calibration were “very good”. The value of NSE for validation (0.87) was also in the “very good” range. Figure 8b shows that some of the predicted groundwater depths were underestimated, whereas others were overestimated. A mean bias error (MBE) value of 0.89 m indicates that, as a whole, the groundwater levels are slightly overestimated. So, one has to be careful in the judgement of simulations. A very good NSE value does not guarantee perfect simulations. There should at least be one other check (Figure 8b) to decide the fate of the model results.

3.2. ANFIS Results

Figure 8a,b presents the results of ANFIS. It is observed that ANFIS is highly efficient technique for simulation of groundwater levels. In Figure 8a,b, the NSE is 0.95 and the 45-degree line is very close to the predicted values of groundwater levels by ANFIS, which shows an excellent performance of ANFIS. Our results are in line with the past studies by Emamgholizadeh et al. [32] and Das [33].

3.3. ANN Model Results

In the case of ANN models, comparatively better performance was observed for ANNs with 10 hidden neurons in both the two-layer and three-layer architecture, as compared with those having 5 or 15 neurons (see Figure 9a,b). The value of NSE was comparatively higher when the models were run with 10 neurons. Figure 10a,b illustrates that the ANN performed better in the case with three layers, as compared with the two-layer architecture. For training, testing, and validation, NSE values were found to be 0.97, 0.94, and 0.92, respectively, for the three-layer architecture and 0.96, 0.86, and 0.91, respectively, for the two-layer architecture. The minimum value of NSE in these results was 0.86, which also represents a “very good” performance (NSE range 0.75 to 1 represents “very good” performance [54].

Further, the performance of the Scaled Conjugate Gradient training function was better than that of the Levenberg–Marquardt and Bayesian Regularization training functions for the three-layer architecture with 10 neurons (Figure 11a–d). Figure 11a shows that the predicted groundwater levels from scaled the conjugate gradient training function were comparatively closer to the 45-degree line (the 45-degree line shows the perfect matching between the predicted and measured groundwater-levels). However, as shown by Figure 11b, the results from the other two training functions, namely Levenberg–Marquardt and Bayesian Regularization, were also not bad. Figure 11c shows that the NSE values remained in the range of 0.65 to 0.97. So, the overall performance of all the three models was “good” to “very good”. This is further strengthened by Figure 11d that shows a very low mean square error in the range of 0.02–0.03. Similar results were observed by previous studies [30,41,42,49]. This is most probably because a stronger optimization scheme was used in the scaled conjugate gradient to get a better global minimum with the minimum possible number of trials [16,54].

3.4. Comparison of ANN, ANFIS and Hydraulic Model Results

Figure 12 presents a comparison of results from various models with respect to their performance. It is observed that the performance of ANN and ANFIS models was better in the training stage as compared to that of the calibration of the MODFLOW. The mean-square-error was lower and the NSE was higher in the training stage of ANN and ANFIS models as compared with the mean-square-error and the NSE in the case of calibration of the MODFLOW. Our results support the results found by recent studies [8,10,32,33,39,55].

3.5. Long-Term Predictions

The predicted groundwater depths for various pumping rates scenarios are shown in Figure 13. The forecast of the drawdown up to the year 2070 with respect to the water level in 2018 was found to be 70–105 m (70, 81.65, 93.3, and 105 m) for four different options of pumping rates (keeping pumping rates constant or increasing by 1%, 2%, and 3% per year as compared to the pumping rate in 2018. The Saq Aquifer is covered by the Qassim Aquifer, which has contaminated shallow groundwater. Such a large drawdown in the Saq aquifer may result in a low pressure in the aquifer, allowing contaminated water from the overlying aquifer to enter into the aquifer. Hence, a serious environmental disaster may occur. So, it is very important to adopt some sustainable management strategies.

This research will be helpful in finding the sustainable pumping rates for future. For example, in case of constant pumping rate, there may be significant increase in the life of the Saq Aquifer. A 33 % reduction in drawdown can be achieved in long-term (50 years in future) by adopting water conservation strategies in comparison to the continuously increasing pumping rates that are required to fulfil the water demands of rapidly increasing population. Results of study will facilitate in the planning and management of limited water resources in Saudi Arabia, the Gulf Region, and other arid regions with similar water stress conditions.

4. Conclusions and Recommendations

The hydraulic model MODFLOW predicted the changes in the groundwater levels of the Saq Aquifer in the vicinity of water supply pumping stations with 55 wells. ANN models effectively simulated the groundwater level changes around the same pumping stations. Two-layered and three-layered ANN models with 5, 10, and 15 hidden neurons in each layer were examined. Three types of training functions were investigated in ANN models. ANFIS has been applied successfully to simulate the aquifer depletion.

The study found that ANFIS can simulate groundwater levels efficiently with NSE up to 0.95. In the case of the ANNs, it provided the best performance with 10 hidden neurons in each layer. The ANNs with 10 hidden neurons in each layer perform better than those with five or 15 neurons. A three-layered architecture with 10 hidden neurons performs better than a two-layered architecture. Moreover, the scaled conjugate gradient performs better than the Levenberg–Marquardt and Bayesian regularization training functions. The hydraulic model is comparatively more reliable, but its performance for the given set of data is not as good as that of the ANN models.

Future predictions show that a substantial groundwater drawdown (ranging from 70 to 150 m) in Saq Aquifer is expected in the next 50 years, with respect to the reference water levels in 2018. There may be a serious environmental disaster if precautionary measures are not taken. However, no increase in pumping rates may enhance the life of the Saq Aquifer significantly. The drawdown will reduce by 33% compared with the drawdown in 2070 in the case of continued pumping at increasing rates to meet the water demands of the increasing population. There is hardly any recharge to the Saq Aquifer and it may vanish after some time. Therefore, new sources of water should be developed as per the goals of the country’s Vision 2030.

Future research work on how new sources of water can be developed is recommended. Improvement in the presentation of the model is also suggested. More effective technology, sophisticated data, knowledge of outcrops, and appropriate management are needed for the study area. Meteorological, hydraulic, hydrologic, and detailed geological data for the whole of the Qassim Region should be collected for wide-ranging hydraulic and ANN modeling. A comparison of various conventional and deep learning artificial intelligence methods is also recommended.