# Simulation of Quantity and Quality of Saq Aquifer Using Artificial Intelligence and Hydraulic Models

Authors: Abdul Razzaq Ghumman; Ghufran Ahmed Pasha; Md. Shafiquzzaman; Afaq Ahmad; Afzal Ahmed; Riaz Akhtar Khan; Rashid Farooq
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Scarcity of water resources is becoming a threatening issue in arid regions like Gulf. Accurate prediction of quantities and quality of groundwater is the first step towards better management of water resources where groundwater is the major source of water supply. Groundwater modelling with respect to its quantity and quality has been performed in this paper using Artificial Neural Networks (ANNs), Adaptive Neurofuzzy Inference System (ANFIS), and hydraulic model MODFLOW. Five types of ANN models with various training functions have been investigated to find the most efficient training function for the prediction of quantity and quality of groundwater, which is an original contribution useful for engineering sector. The results of the hydraulic model, ANFIS, and ANN have been compared. Nash-Sutcliffe Model Efficiency and Mean Square Error have been used for assessing the performance of models. Taylor’s Diagram has also been used to compare various models. The part of Saq Aquifer lying in the Qassim Region has been investigated as the study area. Modern tools, including Geographical Information System (GIS) and Digital Elevation Model (DEM) are applied to process the required data for modelling. Climatic, geographical, and quality of groundwater (contaminants) data are obtained from the Ministry of Environment, Water, and Agriculture, Jeddah/Riyadh. ANFIS model is found to be the most efficient for modelling both the quality and quantity of the aquifer. Sensitivity analysis was performed, and then various future scenarios were investigated for sustainable groundwater pumping. The results of the research will be useful for the community and experts working in the field of water resources engineering, planning, and management in arid regions.


---

## Body

## 1. Introduction

Groundwater resources of the Gulf region are under high stress. These are being depleted at an alarming rate of about 6.8 mm/year in some transboundary aquifers like Saq shared by Iraq, Jordan, and Saudi Arabia. This rapid depletion of groundwater resources (Figure1) is due to anthropogenic and climatic factors. The rainfall in the region is small, but the extraction rate of groundwater for irrigation and domestic purposes is very high due to phenomenal increase in the population [1, 2]. For sustainable management and development of water resources around the globe, the future predictions of groundwater levels are of utmost importance [3–6]. The changes in groundwater levels are affected by several factors. The timing, location, and extent of hydrologic responses to human activities or natural events depend on the duration, intensity, and nature of the activities/event affecting groundwater, the subsoil properties, and the interaction of groundwater and surface water. The field surveys or remote sensing data of various aquifer characteristics may provide a theoretical understanding of the groundwater system, but often the measured data is scarce with respect to space and time, and hence the system-understanding remains undefined and inadequate in most cases.Figure 1
Groundwater depletion and pumping rates for Saq Aquifer (Qassim Region).Some awareness of the complex behavior of the groundwater system can be developed through groundwater models. Once a model is developed and designed properly, it can be used to predict the future behavior of groundwater, which may help in decision-making and exploring the alternative approaches for management, planning, and development of water resources [6–9]. This paper has multiple goals. One of the main goals of this paper is to develop groundwater models for the simulation of water levels of the Saq Aquifer. However, there is no single model which may provide results to the entire satisfaction of the researcher because there is always some degree of uncertainty in model results. Hence a variety of data-driven models based on Artificial Neural Networks (ANNs), Adaptive Neurofuzzy Inference System (ANFIS), and hydraulic model MODFLOW have been investigated in this research.The hydraulic models for transboundary aquifers having undefined boundary conditions for specific regions demand very expensive data acquisition schemes. Although this difficulty has been resolved in this paper by using the “general head boundary condition,” there is a need to investigate some other alternative modelling techniques. The data-driven models can be a good alternative for predicting groundwater levels in the vicinity of some crucial areas like pumping stations [10, 11]. Both types of models, that is, hydraulic and artificial intelligence, are compared in Table 1. Data-driven models are divided into several categories. Khedri et al. [17] compared five different data-driven models for the prediction of groundwater levels. They have investigated the performance of ANN, Adaptive Neurofuzzy Inference System, Fuzzy Logic, Support Vector Regression, and Group Method of Data Handling. According to them, the ANN are the most commonly used efficient method. Majumder and Eldho [18] applied ANN coupled with an optimization technique to model groundwater flow. They have also observed better convergence by ANN combined with optimization. Suprayogi et al. [19] applied ANN to forecast groundwater levels in Meskom, Indonesia. They have observed high values of correlation coefficient having a successful application of ANN for groundwater level forecast for a short time basis. Wunsch et al. [20] have outlined the ability of prediction of groundwater levels by recurrent ANNs using data from t 17 groundwater wells within the Upper Rhine Graben, one of the largest groundwater resources in Central Europe. Solgi et al. [21] have investigated the performance of ANN for long- and short-term groundwater level predictions using more than 17 years of the most recent observed data of groundwater levels in San Antonio, Texas, United States of America. In addition to the prediction of groundwater levels, another dimension of groundwater modelling is predicting its quality [9, 22, 23]. Fadipe et al. [24] have studied groundwater quality for the community of Boluwaduro, Ofatedo, Osun State. The concentrations of nitrite, nitrate, iron, and lead have been investigated by training ANN structure to get the best possible output. Al-Adhaileh and Alsaade [25] have studied groundwater quality for various locations of aquifers in India using various types of artificial intelligence methods. A thorough literature survey shows that there are several dimensions of data-driven models which need to be investigated to choose the best possible model for groundwater predictions.Table 1
Review of data-driven models and hydraulic model MODFLOW.
CharacteristicsANN/ANFISHydraulic modelRemarksEfforts/costs involved in computation expenses (need high-speed computers and time for computations)xAbout 10 timesxAlthough the hydraulic models like MODFLOW require computers with relatively large processing speed as compared to the data-driven models like ANN and ANFIS, which require only ordinary-type computers [12], this factor is not considered nowadays as high-speed computers are easily available.Requirements related to the processing of models and dataxAbout 3 timesxData for only water levels, pumping rates, and concentration of contaminants are required. The data for hydraulic parameters and topography (Digital Elevation Model (DEM)) are required in addition to data for water levels, pumping rates, and concentration of contaminants. The hydraulic model (MODFLOW) requires highly accurate aquifer-parameter values. Furthermore, a mesh of an extraordinary resolution and lower time steps but qualifying the limiting ratio of time step to the nodal distance is required, whereas only a good dataset can be sufficient for data-driven models like ANN and ANFIS [13–15].Model categoryData-driven (black box/semiblack box)Distributed (based on laws of physics)Hydraulic models incorporate physical processes and equations based on laws of physics in predicting groundwater levels and groundwater contamination parameters, whereas ANFIS/ANN, being data-driven models, do not use equations based on laws of physics. These models only use the recorded data for their training, testing, and validation [13–16].Model prediction biasesxAbout 10 timesxThere is a wide range of variation in biases in the case of hydraulic model results, whereas the biases are mostly limited to a certain range for data-driven model (ANFIS/ANN) predictions. However, the ANFIS/ANN lack generality [13–16].Future predictions for a long timeHighly challengingMuch easier once the model is calibratedAfter calibration and validation of the hydraulic model, it is very easy to use it for future predictions, whereas very high experience and expertise in definite phenomenon is required for long-time future predictions by ANFIS/ANN [13–16].In this paper, the ANN techniques and Adaptive Neurofuzzy Inference System (ANFIS) have been investigated rigorously for the prediction of both the quantity and quality of flow through an aquifer. The areas of research regarding ANN include the selection of architecture of ANN, the study of the impact of activation functions, and the choice of an efficient training algorithm. The ANN may have several types of architectures depending upon the number of hidden layers, the number of neurons in each hidden layer, and the activation function [2, 26]. Working of these models further changes if different types of training functions are used [2, 20]. There are twelve types of training functions that can be used. The performances of three types of ANN models in order to predict groundwater levels have been investigated by Coulibaly et al. [27]. Input delay neural network, recurrent neural network, and generalized radial basis function models have been tested and it was found that recurrent neural networks are the most efficient models for predicting monthly groundwater level fluctuations in the Gondo plain, northwestern Burkina Faso. Two types of architectures of ANN have been applied to study groundwater modelling by Nordin et al. [28]. Quantities of different effluents in groundwater have been predicted. The model was applied to real data from groundwater in Faisalabad, the largest industrial city of Pakistan. Seven categories of ANN architectures and training algorithms have been applied by Daliakopoulos et al. [29] to forecast monthly groundwater in Messara Valley in Crete, Greece. They observed that the standard feedforward ANN with the Levenberg-Marquardt algorithm provided the best results for forecasts up to 18 months. Yoon et al. [30] applied ANN with the backpropagation algorithm as a training function to forecast groundwater levels in the beach of the coastal town of Mukho, Donghae city, Korea, on a short-term basis. They found some degree of uncertainty in ANN results and suggested that these should always be designed very carefully. Lohani and Krishan [31] also used backpropagation algorithm as training function in ANN to simulate groundwater level changes in Punjab, India districts of Patiala, Faridkot, Ludhiana, and Ferozepur. They found that groundwater levels can be successfully predicted with acceptable accuracy for future six months using ANN-backpropagation algorithm as training function. The process of learning in advanced types of ANN is achieved by the use of a training algorithm which is similar to optimization [22]. The main aim in each iteration of ANN learning/training is to minimize the error between the predicted and observed groundwater levels and concentration of contaminants.The weights and biases are adjusted in different trials of training in an efficient way to obtain a global minimum with the help of the training algorithm. There are about 12 different training functions for the multilayer perceptron. The ANN will achieve the best performance in the minimum possible computer time if the most efficient training function is used. Performances of only a few out of 12 have been reported in past studies. Hence, in this paper, we have compared the performances of the five most important training functions to simulate the levels and quality of groundwater around pumping stations of Buraydah, Qassim, Saudi Arabia.There are several parameters/contaminants which need to be investigated for studying the quality of the groundwater. However, the main objective of this paper is to compare the performances of various ANN and ANFIS models, so only a few quality parameters have been chosen to make the comparison visible and comparatively easier. TDS, Fe, and Mn were easily available from recorded data. Therefore, these three contaminants have been selected for the quality study. Another point is worth mentioning here that the available data consists of annual values of concentration and groundwater levels. Hence, the investigation regarding seasonal variability is out of the scope of this paper.


## 2. Materials and Methods

### 2.1. Methodological Framework

In this paper, two types of techniques, the hydraulic model MODFLOW and the data-driven models (ANFIS and ANN), have been used to predict the groundwater quality and quantity for the Saq Aquifer. A conceptual comparison of the two techniques is given in Table1. An overall methodology framework is shown in Figure 2. Three main steps are involved in this research work. The first step comprised data collection, including pumping rates, groundwater levels, well locations, the concentration of contaminants, and parameters of the aquifer. The data for pumping rates, groundwater levels, and concentration of contaminants have been normalized in the case of ANFIS and ANN models. The second step was to calibrate/train/test/validate as per the requirement of both categories of models (hydraulic/ANFIS andANN). Part of the data (say 60%) was used for calibration/training, and the remaining (40%) was used for validation/testing. It is worth mentioning that calibration and validation were done for both the quantity and quality of groundwater. Hydraulic parameters were adjusted to obtain accurate groundwater levels at certain locations after a certain model runtime as close to the observed groundwater levels as possible. The calibration of the pollutant transport model has mainly been done by adjusting the contaminant transport parameters (dispersion coefficient, etc.) until the observed concentration values were consistent with the model's predicted values. Six wells were processed and used as concentration observation wells with their coordinates and historical records of groundwater levels.Figure 2
Groundwater modelling framework to explore the impacts of pumping on depletion of aquifer.The data were divided into two parts, 60% for calibration and 40% for validation of the hydraulic model, and three parts, 60%, 20%, and 20%, for the training, validation, and testing of ANFIS/ANN models, respectively. A computer program (coding) was prepared in MATLAB to develop ANFIS/ANN models to achieve the required objectives.Sensitivity analysis of parameters was made to see the impact of various hydraulic and quality parameters on model results. Step three included the application of the chosen model for future predictions of groundwater quality and levels to study the aquifer depletion in both quality and quantity. Due to the phenomenal increase in population and high living standards of people, there is an increasing trend in groundwater use [32]. Eight future scenarios for different pumping rates (as given below) have been tested to investigate the depletion of an aquifer with respect to both quality and quantity. The farthest optimistic scenario, farthest pessimistic scenario, and possible scenarios between the two extremes have been examined as given in Table 2.Table 2
Future scenarios for pumping rates.
Scenario no.DescriptionScenario no.Description1Assuming a decrease in the pumping rate by 1% yearly5Assuming the pumping rate to increase 3% yearly2Assuming the present pumping rate to continue during the entire period6Assuming the pumping rate to increase 4% yearly3Assuming the pumping rate to increase 1% yearly7Assuming the pumping rate to increase 5% yearly4Assuming the pumping rate to increase 2% yearly8Assuming scenario 2 with an artificial recharge case


### 2.2. Study Area

The study area included the most important area of Saq Aquifer located in the Qassim Region, as shown in Figure3. Qassim region is one of the important regions of the Kingdom of Saudi Arabia. It is characterized as an arid region. It has an area of approximately 80,000 km2, located between latitudes 25°00′ and 27°00′ N and longitudes 42°30′ and 45°00′ E [2, 5]. It has an altitude ranging between 600 and 850 m above sea level. Qassim has a population of nearly 1,423,935 capita as per 2017 records. It is the second-largest agricultural region in Saudi Arabia, and groundwater from the Saq Aquifer is available easily. The region comprises a desert climate, characterized by cold and rainy winter and hot summer with low humidity. The evaporation is nearly 4500 mm per annum on average. The summer temperature is as high as 49°C during the daytime and 36°C during the nighttime. However, the winter temperature may sometimes fall below 0°C. The average annual rainfall in Qassim Region is about 125 mm [2]. The main source of water supply to the region is groundwater extracted mainly from the Saq Aquifer, which is one of the most important aquifers in the Kingdom [5, 33]. It has a very vast outcrop, which spreads over a length of approximately 1200 km only in the Kingdom of Saudi Arabia and joins the border of Jordan as far south as latitude 24°30′ N and longitude 45°00′ E. The huge aquifer is confined in nature with a thickness varying from 400 m in the southern part to 700 m in its northern part. The part of the aquifer under study has an average thickness of about 500 m [2]. The strata consist of sandstone having medium to course size. A few areas comprise sandstone of fine size at local levels. Most of the private and public tube wells are partially penetrating, with 125 m screen and a depth of approximately 650 m [2]. Wadi Al-Rummah is considered a big valley that passes through the Saq area in the Qassim Region. The Saq Aquifer has a high salinity value in the area near Wadi Al-Rummah compared to the area away from Wadi Al-Rummah. Recharge from Wadi, agriculture processes, and effluent from wastewater treatment plants are essential factors that may contribute to its comparatively higher salinity in the surrounding area of Wadi Al-Rummah [5].Figure 3
Map of the study area showing the locations of pumping wells, Saq in the Qassim Region, and Saq outcrop.


### 2.3. Description of Hydraulic Model (MODFLOW)

#### 2.3.1. Groundwater Levels and Flow Modelling

The groundwater flow modelling is highly complicated [34] as it is governed by partial differential equations as given below [35].(1)∂∂xKxx∂h∂x+∂∂yKyy∂h∂y+∂∂zKzz∂h∂z=Ss∂h∂t−W,where the parameters “Kxx,” “Kyy,” and “Kzz” represent the hydraulic conductivities in the x, y, and z directions, respectively; the groundwater head is expressed by “h”; “Ss” is the parameter called the specific storage of the aquifer; the source/sink is represented by W; and “t” is time. A denser mesh was created at the locations of wells to model the wells as accurately as possible [36–39].Visual MODFLOW is a commercial graphical user interface. It has unique features in which users have the ability to use Excel files, Surfer grids, GIS, and AutoCAD data as input files instead of the form of text files. Also, the results can be shown in the charts and contour maps. These features help in lowering the execution time and making the analysis of results easier [37–39]. In addition to finite element, there are many finite difference schemes to solve (1), including explicit finite difference and implicit finite difference. The model MODFLOW uses an implicit finite difference scheme to achieve the solution of equation (1). The input and output data for MODFLOW are presented in Table 3. Selection of the mesh size and appropriate boundary conditions is a highly challenging task in large aquifers like Saq. In this paper, a rectangular type of mesh of 30 km × 20 km was developed by creating 80 rows and 70 columns; in one case, slightly bigger area was modelled in case of studying the area around Wadi Al-Rummah. The numbers of columns and rows for the mesh around the pumping wells were significantly increased to determine a highly dense mesh for the prediction of the quantity and quality of groundwater as accurately as possible. The thickness of the aquifer in the vertical direction was taken as approximately 500 m with an impermeable cover of about 628 m above 500 m thick layer of Saq Aquifer [2, 33]. As mentioned above, the pumping wells are partially penetrating with a 125 m long screen. Therefore, for simulating the effects of well screen accurately, the confined 500 m thick aquifer has been divided into four layers of equal thickness (125 m each). Accordingly, five layers of 631, 125, 125, 125, and 125 m thickness have been adopted (the top 631 m layer being impermeable) (Figure 4).Table 3
Input values of data for MODFLOW.
Required dataData description/estimationAdopted values of dataTransmissivity and conductivityCooper-Jacob technique applied on pumping test to estimate the hydraulic conductivity of Saq Aquifer in Qassim Region [33, 40]T = 0.024 to 1.62 m2/minute, Kyy = Kxx = Kzz = 2.8 m/dayNo. of aquifer layersPublished research papers [2, 33]Adopted layers = 5, each of 631, 125, 125, 125, and 125 m, respectivelyMesh sizeBy running model with different mesh sizes and choosing the bestAs much dense as possible. It was further made dense around the wellsPumping ratesAs per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaPublic and private tube wells data available from 1980 to 2018Water levels (a)As per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaPublic and private tube wells data available from 1980 to 2018Initial conditionsAs per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaWater levels at the start time of simulationsBoundary conditionsCalibrated assuming general head boundary conditionsHighly challenging workStorage coefficientCooper-Jacob technique applied on pumping test to estimate the hydraulic conductivity of Saq Aquifer in Qassim Region [33, 40]2.5 × 10 – 3 to 6.3 × 10 – 5Contaminant concentrationsRecorded valuesRecorded dataFigure 4
Layers of the Saq Aquifer for MODFLOW models (after Almuhaylan et al. [2]).


#### 2.3.2. Quality Modelling

For groundwater quality investigations, the transport and concentration of contaminants are modelled. Hence an additional equation is required for quality modelling on the basis of the transport of contaminants. The equation given below is a partial differential equation for the three-dimensional transport of contaminants in a groundwater flow system [41]:(2)∂∁∂t=∂∂xiDij∂∁∂xi−∂∂xiυiC+qsθCs+∑1NRk,where “C” is the dissolved concentration of contaminant k in groundwater, “θ” is porosity of the subsurface porous medium, “Dij” is hydrodynamic dispersion coefficient, “xi” is the distance along the respective Cartesian coordinate axis, “vi” is seepage or linear pore water velocity, “qs” is volumetric flow rate, “Rk” is the chemical reaction term, and “Cs” is the concentration of the source [41, 42].There are several techniques to solve (2) combined with equation (1). Finite difference and finite element are the most commonly used techniques. Simultaneous solutions of equations (1) and (2) by numerical methods are used for groundwater quantity and quality modelling, which is a challenging task. Several commercial packages are available for this purpose. MODFLOW is used in this paper to model both the quality and quantity of groundwater.All the above parameters described under Section2.3.1 and used in the prediction of groundwater levels were kept the same for contaminant transport simulations (quality modelling). The constant concentration option has been selected as a transport boundary condition in the Visual MODFLOW. Average values of TDS, Mn, and Fe during 2010 were taken as initial values in the model. To assign initial concentration values, “Assign” from the menu bar of MODFLOW was selected, and the values of initial concentration were added in “Assign Window.” The MT3DMS module of Visual MODFLOW was used to simulate the contaminants transport route and concentration of contaminants in groundwater. This module solves (2) coupled with (1) by finite difference technique. The standard dispersion parameters have been adopted. The option of saturated constant density was chosen. The calibration of the contaminant (TDS, Fe, and Mn) transport model was executed by adjusting the contaminant transport parameters (dispersion coefficient, etc.) until the observed concentration values of TDS, Fe, and Mn were consistent with the predicted values of the model for these contaminants. As described above, the six wells were used as the observation wells for concentration. The coordinates of these wells, concentration of TDS, Fe, and Mn, and historical records of groundwater levels were provided to the model.


### 2.4. Artificial Intelligence Models

In groundwater modelling, the application of models based on artificial intelligence (AI) is becoming highly popular and is useful in fact [2, 20, 22, 26]. AI techniques have been effectively used in water resources planning/development/management, hydrological processes, water quality predictions, and operation of reservoirs [43–46]. The AI models are gaining importance and expansion in application scope because of their ability to generalize the relationships for a complex phenomenon based on the learning process [47, 48]. However, the use of AI models in the area of specialization of prediction of both groundwater quality and quantity is limited and is still a vital and critical research issue [47–49]. A general layout of AI models is shown in Figures 5(a) and 5(b).Figure 5
AI modelling approach. (a) Flowchart representing the MLP (a feedforward-ANN). (b) Five functional blocks in ANFIS.
(a)(b)


#### 2.4.1. Artificial Neural Network (ANN) Models

In ANN type models, the neurons play a key role in determining the values of output variables comparable to the target values from input variable values. As shown in Figure5(a), there are mainly three layers: input, hidden, and output layers. An internal-structural procedure is performing the calculations. The hidden layers in an internal-structural procedure, that is, the backbone of the predictions, work as the backpropagation (BP) and the feedforward mechanism for delivering outputs with required precision. The case of multiple hidden layers generates the well-known “Multilayer Perceptron” (MPL). The training function is another important entity of ANN models which are of several types [50]. This paper has used 5 distinct training functions (Table 4). On these bases, the ANN models are termed as ANN1, ANN2, ANN3, ANN4, and ANN5. Both the triple and double hidden layers have been tested. These are represented by the symbol “DL” for the double layer and the symbol “TL” for the triple layer. The hidden layers may have a different number of neurons which need to be fixed for the best possible results [51–53]. The hit-and-trial method has been used in this research for the selection of the best possible number of neurons. A model is considered to be precise if its predictions are matching to the measured values. In modern techniques of ANN, it is accomplished by combining an effective optimization methodology in which weights are adjusted ANN process to find the highest possible accuracy.Table 4
Various training functions used for ANN models.
ModelDescription of training functionFunctionANN1Scaled Conjugate Gradient BPTrainscgANN2Levenberg-Marquardt BPTrainlmANN3Bayesian regularizationTrainbrANN4BFGS Quasi-Newton BPTrainbfgANN5Resilient backpropagation (rprop)Trainrp


#### 2.4.2. Adaptive Neurofuzzy Inference Systems

An effective optimization subroutine with a hybrid AI model may become a basis for a robust technique for modelling the quality and quantity of groundwater. It is termed as an “Adaptive Neurofuzzy Inference System” (ANFIS). ANFIS can be used to solve a complex nonlinear hydraulic phenomenon with high precision [2, 50]. Figure 5(b) shows the structure of ANFIS, including five layers (Table 5).Table 5
Structural layers of ANFIS.
Layer no.Description1A fuzzy layer in which every intersection is called an adaptive node. The member functions may be treated as nodal functions like the generalized bell, and Gaussian functions may be used as nodal functions2A product layer in which every nodal output represents the firing strength of a rule3A normalized layer in which each node expresses the firing strength of a rule in the normalized form4A defuzzy layer in which each node represents an adaptive node with a nodal function specifying the influence of the rules in terms of the overall output5A total output layer that describes the defuzzified output values


### 2.5. Model Performance Evaluation

The performances of models in this paper are represented by two statistical parameters, the Nash-Sutcliffe Model Efficiency (NSE) and Mean Square Error (MSE). Both parameters are based on some weighted difference between the recorded and simulated values of variables (groundwater levels, TDS, Fe, and Mn). About 38 years (1980–2018) of data for groundwater levels and pumping rates and long records of TDS, Mn, and Fe have been used to evaluate the model performance for different runtime stages (training/calibration, testing, and validation). Many of the previous studies have utilized a similar approach [2, 50, 54].The Nash-Sutcliffe Model Efficiency (NSE) is given as(3)NSE=1−∑i=1nGWLio−GWLip2∑i=1nGWLio−GWLavg2,where “GWLio and GWLip)2” denote the observed and predicted groundwater levels, respectively, for the ith data value and “n” is the total number of observed/predicted data points. The NSE values ranging from 0.75 to 1.0 may be termed as “very good,” 0.65 to 0.75 may be taken as “good,” and 0.5 to 0.65 as “satisfactory,” and the values ranging between 0.4 and 0.5 may be expressed as “acceptable” regarding the model’s performance [55, 56].The Mean Square Error (MSE) is given as(4)MSE=1n∑i=1nGWLio−GWLip2.In the case of quality modelling,GWL will be replaced by the values of TDS or the concentrations of Mn and Fe in (3) and (4).


## 2.1. Methodological Framework

In this paper, two types of techniques, the hydraulic model MODFLOW and the data-driven models (ANFIS and ANN), have been used to predict the groundwater quality and quantity for the Saq Aquifer. A conceptual comparison of the two techniques is given in Table1. An overall methodology framework is shown in Figure 2. Three main steps are involved in this research work. The first step comprised data collection, including pumping rates, groundwater levels, well locations, the concentration of contaminants, and parameters of the aquifer. The data for pumping rates, groundwater levels, and concentration of contaminants have been normalized in the case of ANFIS and ANN models. The second step was to calibrate/train/test/validate as per the requirement of both categories of models (hydraulic/ANFIS andANN). Part of the data (say 60%) was used for calibration/training, and the remaining (40%) was used for validation/testing. It is worth mentioning that calibration and validation were done for both the quantity and quality of groundwater. Hydraulic parameters were adjusted to obtain accurate groundwater levels at certain locations after a certain model runtime as close to the observed groundwater levels as possible. The calibration of the pollutant transport model has mainly been done by adjusting the contaminant transport parameters (dispersion coefficient, etc.) until the observed concentration values were consistent with the model's predicted values. Six wells were processed and used as concentration observation wells with their coordinates and historical records of groundwater levels.Figure 2
Groundwater modelling framework to explore the impacts of pumping on depletion of aquifer.The data were divided into two parts, 60% for calibration and 40% for validation of the hydraulic model, and three parts, 60%, 20%, and 20%, for the training, validation, and testing of ANFIS/ANN models, respectively. A computer program (coding) was prepared in MATLAB to develop ANFIS/ANN models to achieve the required objectives.Sensitivity analysis of parameters was made to see the impact of various hydraulic and quality parameters on model results. Step three included the application of the chosen model for future predictions of groundwater quality and levels to study the aquifer depletion in both quality and quantity. Due to the phenomenal increase in population and high living standards of people, there is an increasing trend in groundwater use [32]. Eight future scenarios for different pumping rates (as given below) have been tested to investigate the depletion of an aquifer with respect to both quality and quantity. The farthest optimistic scenario, farthest pessimistic scenario, and possible scenarios between the two extremes have been examined as given in Table 2.Table 2
Future scenarios for pumping rates.
Scenario no.DescriptionScenario no.Description1Assuming a decrease in the pumping rate by 1% yearly5Assuming the pumping rate to increase 3% yearly2Assuming the present pumping rate to continue during the entire period6Assuming the pumping rate to increase 4% yearly3Assuming the pumping rate to increase 1% yearly7Assuming the pumping rate to increase 5% yearly4Assuming the pumping rate to increase 2% yearly8Assuming scenario 2 with an artificial recharge case


## 2.2. Study Area

The study area included the most important area of Saq Aquifer located in the Qassim Region, as shown in Figure3. Qassim region is one of the important regions of the Kingdom of Saudi Arabia. It is characterized as an arid region. It has an area of approximately 80,000 km2, located between latitudes 25°00′ and 27°00′ N and longitudes 42°30′ and 45°00′ E [2, 5]. It has an altitude ranging between 600 and 850 m above sea level. Qassim has a population of nearly 1,423,935 capita as per 2017 records. It is the second-largest agricultural region in Saudi Arabia, and groundwater from the Saq Aquifer is available easily. The region comprises a desert climate, characterized by cold and rainy winter and hot summer with low humidity. The evaporation is nearly 4500 mm per annum on average. The summer temperature is as high as 49°C during the daytime and 36°C during the nighttime. However, the winter temperature may sometimes fall below 0°C. The average annual rainfall in Qassim Region is about 125 mm [2]. The main source of water supply to the region is groundwater extracted mainly from the Saq Aquifer, which is one of the most important aquifers in the Kingdom [5, 33]. It has a very vast outcrop, which spreads over a length of approximately 1200 km only in the Kingdom of Saudi Arabia and joins the border of Jordan as far south as latitude 24°30′ N and longitude 45°00′ E. The huge aquifer is confined in nature with a thickness varying from 400 m in the southern part to 700 m in its northern part. The part of the aquifer under study has an average thickness of about 500 m [2]. The strata consist of sandstone having medium to course size. A few areas comprise sandstone of fine size at local levels. Most of the private and public tube wells are partially penetrating, with 125 m screen and a depth of approximately 650 m [2]. Wadi Al-Rummah is considered a big valley that passes through the Saq area in the Qassim Region. The Saq Aquifer has a high salinity value in the area near Wadi Al-Rummah compared to the area away from Wadi Al-Rummah. Recharge from Wadi, agriculture processes, and effluent from wastewater treatment plants are essential factors that may contribute to its comparatively higher salinity in the surrounding area of Wadi Al-Rummah [5].Figure 3
Map of the study area showing the locations of pumping wells, Saq in the Qassim Region, and Saq outcrop.


## 2.3. Description of Hydraulic Model (MODFLOW)

### 2.3.1. Groundwater Levels and Flow Modelling

The groundwater flow modelling is highly complicated [34] as it is governed by partial differential equations as given below [35].(1)∂∂xKxx∂h∂x+∂∂yKyy∂h∂y+∂∂zKzz∂h∂z=Ss∂h∂t−W,where the parameters “Kxx,” “Kyy,” and “Kzz” represent the hydraulic conductivities in the x, y, and z directions, respectively; the groundwater head is expressed by “h”; “Ss” is the parameter called the specific storage of the aquifer; the source/sink is represented by W; and “t” is time. A denser mesh was created at the locations of wells to model the wells as accurately as possible [36–39].Visual MODFLOW is a commercial graphical user interface. It has unique features in which users have the ability to use Excel files, Surfer grids, GIS, and AutoCAD data as input files instead of the form of text files. Also, the results can be shown in the charts and contour maps. These features help in lowering the execution time and making the analysis of results easier [37–39]. In addition to finite element, there are many finite difference schemes to solve (1), including explicit finite difference and implicit finite difference. The model MODFLOW uses an implicit finite difference scheme to achieve the solution of equation (1). The input and output data for MODFLOW are presented in Table 3. Selection of the mesh size and appropriate boundary conditions is a highly challenging task in large aquifers like Saq. In this paper, a rectangular type of mesh of 30 km × 20 km was developed by creating 80 rows and 70 columns; in one case, slightly bigger area was modelled in case of studying the area around Wadi Al-Rummah. The numbers of columns and rows for the mesh around the pumping wells were significantly increased to determine a highly dense mesh for the prediction of the quantity and quality of groundwater as accurately as possible. The thickness of the aquifer in the vertical direction was taken as approximately 500 m with an impermeable cover of about 628 m above 500 m thick layer of Saq Aquifer [2, 33]. As mentioned above, the pumping wells are partially penetrating with a 125 m long screen. Therefore, for simulating the effects of well screen accurately, the confined 500 m thick aquifer has been divided into four layers of equal thickness (125 m each). Accordingly, five layers of 631, 125, 125, 125, and 125 m thickness have been adopted (the top 631 m layer being impermeable) (Figure 4).Table 3
Input values of data for MODFLOW.
Required dataData description/estimationAdopted values of dataTransmissivity and conductivityCooper-Jacob technique applied on pumping test to estimate the hydraulic conductivity of Saq Aquifer in Qassim Region [33, 40]T = 0.024 to 1.62 m2/minute, Kyy = Kxx = Kzz = 2.8 m/dayNo. of aquifer layersPublished research papers [2, 33]Adopted layers = 5, each of 631, 125, 125, 125, and 125 m, respectivelyMesh sizeBy running model with different mesh sizes and choosing the bestAs much dense as possible. It was further made dense around the wellsPumping ratesAs per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaPublic and private tube wells data available from 1980 to 2018Water levels (a)As per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaPublic and private tube wells data available from 1980 to 2018Initial conditionsAs per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaWater levels at the start time of simulationsBoundary conditionsCalibrated assuming general head boundary conditionsHighly challenging workStorage coefficientCooper-Jacob technique applied on pumping test to estimate the hydraulic conductivity of Saq Aquifer in Qassim Region [33, 40]2.5 × 10 – 3 to 6.3 × 10 – 5Contaminant concentrationsRecorded valuesRecorded dataFigure 4
Layers of the Saq Aquifer for MODFLOW models (after Almuhaylan et al. [2]).


### 2.3.2. Quality Modelling

For groundwater quality investigations, the transport and concentration of contaminants are modelled. Hence an additional equation is required for quality modelling on the basis of the transport of contaminants. The equation given below is a partial differential equation for the three-dimensional transport of contaminants in a groundwater flow system [41]:(2)∂∁∂t=∂∂xiDij∂∁∂xi−∂∂xiυiC+qsθCs+∑1NRk,where “C” is the dissolved concentration of contaminant k in groundwater, “θ” is porosity of the subsurface porous medium, “Dij” is hydrodynamic dispersion coefficient, “xi” is the distance along the respective Cartesian coordinate axis, “vi” is seepage or linear pore water velocity, “qs” is volumetric flow rate, “Rk” is the chemical reaction term, and “Cs” is the concentration of the source [41, 42].There are several techniques to solve (2) combined with equation (1). Finite difference and finite element are the most commonly used techniques. Simultaneous solutions of equations (1) and (2) by numerical methods are used for groundwater quantity and quality modelling, which is a challenging task. Several commercial packages are available for this purpose. MODFLOW is used in this paper to model both the quality and quantity of groundwater.All the above parameters described under Section2.3.1 and used in the prediction of groundwater levels were kept the same for contaminant transport simulations (quality modelling). The constant concentration option has been selected as a transport boundary condition in the Visual MODFLOW. Average values of TDS, Mn, and Fe during 2010 were taken as initial values in the model. To assign initial concentration values, “Assign” from the menu bar of MODFLOW was selected, and the values of initial concentration were added in “Assign Window.” The MT3DMS module of Visual MODFLOW was used to simulate the contaminants transport route and concentration of contaminants in groundwater. This module solves (2) coupled with (1) by finite difference technique. The standard dispersion parameters have been adopted. The option of saturated constant density was chosen. The calibration of the contaminant (TDS, Fe, and Mn) transport model was executed by adjusting the contaminant transport parameters (dispersion coefficient, etc.) until the observed concentration values of TDS, Fe, and Mn were consistent with the predicted values of the model for these contaminants. As described above, the six wells were used as the observation wells for concentration. The coordinates of these wells, concentration of TDS, Fe, and Mn, and historical records of groundwater levels were provided to the model.


## 2.3.1. Groundwater Levels and Flow Modelling

The groundwater flow modelling is highly complicated [34] as it is governed by partial differential equations as given below [35].(1)∂∂xKxx∂h∂x+∂∂yKyy∂h∂y+∂∂zKzz∂h∂z=Ss∂h∂t−W,where the parameters “Kxx,” “Kyy,” and “Kzz” represent the hydraulic conductivities in the x, y, and z directions, respectively; the groundwater head is expressed by “h”; “Ss” is the parameter called the specific storage of the aquifer; the source/sink is represented by W; and “t” is time. A denser mesh was created at the locations of wells to model the wells as accurately as possible [36–39].Visual MODFLOW is a commercial graphical user interface. It has unique features in which users have the ability to use Excel files, Surfer grids, GIS, and AutoCAD data as input files instead of the form of text files. Also, the results can be shown in the charts and contour maps. These features help in lowering the execution time and making the analysis of results easier [37–39]. In addition to finite element, there are many finite difference schemes to solve (1), including explicit finite difference and implicit finite difference. The model MODFLOW uses an implicit finite difference scheme to achieve the solution of equation (1). The input and output data for MODFLOW are presented in Table 3. Selection of the mesh size and appropriate boundary conditions is a highly challenging task in large aquifers like Saq. In this paper, a rectangular type of mesh of 30 km × 20 km was developed by creating 80 rows and 70 columns; in one case, slightly bigger area was modelled in case of studying the area around Wadi Al-Rummah. The numbers of columns and rows for the mesh around the pumping wells were significantly increased to determine a highly dense mesh for the prediction of the quantity and quality of groundwater as accurately as possible. The thickness of the aquifer in the vertical direction was taken as approximately 500 m with an impermeable cover of about 628 m above 500 m thick layer of Saq Aquifer [2, 33]. As mentioned above, the pumping wells are partially penetrating with a 125 m long screen. Therefore, for simulating the effects of well screen accurately, the confined 500 m thick aquifer has been divided into four layers of equal thickness (125 m each). Accordingly, five layers of 631, 125, 125, 125, and 125 m thickness have been adopted (the top 631 m layer being impermeable) (Figure 4).Table 3
Input values of data for MODFLOW.
Required dataData description/estimationAdopted values of dataTransmissivity and conductivityCooper-Jacob technique applied on pumping test to estimate the hydraulic conductivity of Saq Aquifer in Qassim Region [33, 40]T = 0.024 to 1.62 m2/minute, Kyy = Kxx = Kzz = 2.8 m/dayNo. of aquifer layersPublished research papers [2, 33]Adopted layers = 5, each of 631, 125, 125, 125, and 125 m, respectivelyMesh sizeBy running model with different mesh sizes and choosing the bestAs much dense as possible. It was further made dense around the wellsPumping ratesAs per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaPublic and private tube wells data available from 1980 to 2018Water levels (a)As per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaPublic and private tube wells data available from 1980 to 2018Initial conditionsAs per data collected from the Ministry of Environment, Water, and Agriculture, Saudi ArabiaWater levels at the start time of simulationsBoundary conditionsCalibrated assuming general head boundary conditionsHighly challenging workStorage coefficientCooper-Jacob technique applied on pumping test to estimate the hydraulic conductivity of Saq Aquifer in Qassim Region [33, 40]2.5 × 10 – 3 to 6.3 × 10 – 5Contaminant concentrationsRecorded valuesRecorded dataFigure 4
Layers of the Saq Aquifer for MODFLOW models (after Almuhaylan et al. [2]).


## 2.3.2. Quality Modelling

For groundwater quality investigations, the transport and concentration of contaminants are modelled. Hence an additional equation is required for quality modelling on the basis of the transport of contaminants. The equation given below is a partial differential equation for the three-dimensional transport of contaminants in a groundwater flow system [41]:(2)∂∁∂t=∂∂xiDij∂∁∂xi−∂∂xiυiC+qsθCs+∑1NRk,where “C” is the dissolved concentration of contaminant k in groundwater, “θ” is porosity of the subsurface porous medium, “Dij” is hydrodynamic dispersion coefficient, “xi” is the distance along the respective Cartesian coordinate axis, “vi” is seepage or linear pore water velocity, “qs” is volumetric flow rate, “Rk” is the chemical reaction term, and “Cs” is the concentration of the source [41, 42].There are several techniques to solve (2) combined with equation (1). Finite difference and finite element are the most commonly used techniques. Simultaneous solutions of equations (1) and (2) by numerical methods are used for groundwater quantity and quality modelling, which is a challenging task. Several commercial packages are available for this purpose. MODFLOW is used in this paper to model both the quality and quantity of groundwater.All the above parameters described under Section2.3.1 and used in the prediction of groundwater levels were kept the same for contaminant transport simulations (quality modelling). The constant concentration option has been selected as a transport boundary condition in the Visual MODFLOW. Average values of TDS, Mn, and Fe during 2010 were taken as initial values in the model. To assign initial concentration values, “Assign” from the menu bar of MODFLOW was selected, and the values of initial concentration were added in “Assign Window.” The MT3DMS module of Visual MODFLOW was used to simulate the contaminants transport route and concentration of contaminants in groundwater. This module solves (2) coupled with (1) by finite difference technique. The standard dispersion parameters have been adopted. The option of saturated constant density was chosen. The calibration of the contaminant (TDS, Fe, and Mn) transport model was executed by adjusting the contaminant transport parameters (dispersion coefficient, etc.) until the observed concentration values of TDS, Fe, and Mn were consistent with the predicted values of the model for these contaminants. As described above, the six wells were used as the observation wells for concentration. The coordinates of these wells, concentration of TDS, Fe, and Mn, and historical records of groundwater levels were provided to the model.


## 2.4. Artificial Intelligence Models

In groundwater modelling, the application of models based on artificial intelligence (AI) is becoming highly popular and is useful in fact [2, 20, 22, 26]. AI techniques have been effectively used in water resources planning/development/management, hydrological processes, water quality predictions, and operation of reservoirs [43–46]. The AI models are gaining importance and expansion in application scope because of their ability to generalize the relationships for a complex phenomenon based on the learning process [47, 48]. However, the use of AI models in the area of specialization of prediction of both groundwater quality and quantity is limited and is still a vital and critical research issue [47–49]. A general layout of AI models is shown in Figures 5(a) and 5(b).Figure 5
AI modelling approach. (a) Flowchart representing the MLP (a feedforward-ANN). (b) Five functional blocks in ANFIS.
(a)(b)


### 2.4.1. Artificial Neural Network (ANN) Models

In ANN type models, the neurons play a key role in determining the values of output variables comparable to the target values from input variable values. As shown in Figure5(a), there are mainly three layers: input, hidden, and output layers. An internal-structural procedure is performing the calculations. The hidden layers in an internal-structural procedure, that is, the backbone of the predictions, work as the backpropagation (BP) and the feedforward mechanism for delivering outputs with required precision. The case of multiple hidden layers generates the well-known “Multilayer Perceptron” (MPL). The training function is another important entity of ANN models which are of several types [50]. This paper has used 5 distinct training functions (Table 4). On these bases, the ANN models are termed as ANN1, ANN2, ANN3, ANN4, and ANN5. Both the triple and double hidden layers have been tested. These are represented by the symbol “DL” for the double layer and the symbol “TL” for the triple layer. The hidden layers may have a different number of neurons which need to be fixed for the best possible results [51–53]. The hit-and-trial method has been used in this research for the selection of the best possible number of neurons. A model is considered to be precise if its predictions are matching to the measured values. In modern techniques of ANN, it is accomplished by combining an effective optimization methodology in which weights are adjusted ANN process to find the highest possible accuracy.Table 4
Various training functions used for ANN models.
ModelDescription of training functionFunctionANN1Scaled Conjugate Gradient BPTrainscgANN2Levenberg-Marquardt BPTrainlmANN3Bayesian regularizationTrainbrANN4BFGS Quasi-Newton BPTrainbfgANN5Resilient backpropagation (rprop)Trainrp


### 2.4.2. Adaptive Neurofuzzy Inference Systems

An effective optimization subroutine with a hybrid AI model may become a basis for a robust technique for modelling the quality and quantity of groundwater. It is termed as an “Adaptive Neurofuzzy Inference System” (ANFIS). ANFIS can be used to solve a complex nonlinear hydraulic phenomenon with high precision [2, 50]. Figure 5(b) shows the structure of ANFIS, including five layers (Table 5).Table 5
Structural layers of ANFIS.
Layer no.Description1A fuzzy layer in which every intersection is called an adaptive node. The member functions may be treated as nodal functions like the generalized bell, and Gaussian functions may be used as nodal functions2A product layer in which every nodal output represents the firing strength of a rule3A normalized layer in which each node expresses the firing strength of a rule in the normalized form4A defuzzy layer in which each node represents an adaptive node with a nodal function specifying the influence of the rules in terms of the overall output5A total output layer that describes the defuzzified output values


## 2.4.1. Artificial Neural Network (ANN) Models

In ANN type models, the neurons play a key role in determining the values of output variables comparable to the target values from input variable values. As shown in Figure5(a), there are mainly three layers: input, hidden, and output layers. An internal-structural procedure is performing the calculations. The hidden layers in an internal-structural procedure, that is, the backbone of the predictions, work as the backpropagation (BP) and the feedforward mechanism for delivering outputs with required precision. The case of multiple hidden layers generates the well-known “Multilayer Perceptron” (MPL). The training function is another important entity of ANN models which are of several types [50]. This paper has used 5 distinct training functions (Table 4). On these bases, the ANN models are termed as ANN1, ANN2, ANN3, ANN4, and ANN5. Both the triple and double hidden layers have been tested. These are represented by the symbol “DL” for the double layer and the symbol “TL” for the triple layer. The hidden layers may have a different number of neurons which need to be fixed for the best possible results [51–53]. The hit-and-trial method has been used in this research for the selection of the best possible number of neurons. A model is considered to be precise if its predictions are matching to the measured values. In modern techniques of ANN, it is accomplished by combining an effective optimization methodology in which weights are adjusted ANN process to find the highest possible accuracy.Table 4
Various training functions used for ANN models.
ModelDescription of training functionFunctionANN1Scaled Conjugate Gradient BPTrainscgANN2Levenberg-Marquardt BPTrainlmANN3Bayesian regularizationTrainbrANN4BFGS Quasi-Newton BPTrainbfgANN5Resilient backpropagation (rprop)Trainrp


## 2.4.2. Adaptive Neurofuzzy Inference Systems

An effective optimization subroutine with a hybrid AI model may become a basis for a robust technique for modelling the quality and quantity of groundwater. It is termed as an “Adaptive Neurofuzzy Inference System” (ANFIS). ANFIS can be used to solve a complex nonlinear hydraulic phenomenon with high precision [2, 50]. Figure 5(b) shows the structure of ANFIS, including five layers (Table 5).Table 5
Structural layers of ANFIS.
Layer no.Description1A fuzzy layer in which every intersection is called an adaptive node. The member functions may be treated as nodal functions like the generalized bell, and Gaussian functions may be used as nodal functions2A product layer in which every nodal output represents the firing strength of a rule3A normalized layer in which each node expresses the firing strength of a rule in the normalized form4A defuzzy layer in which each node represents an adaptive node with a nodal function specifying the influence of the rules in terms of the overall output5A total output layer that describes the defuzzified output values


## 2.5. Model Performance Evaluation

The performances of models in this paper are represented by two statistical parameters, the Nash-Sutcliffe Model Efficiency (NSE) and Mean Square Error (MSE). Both parameters are based on some weighted difference between the recorded and simulated values of variables (groundwater levels, TDS, Fe, and Mn). About 38 years (1980–2018) of data for groundwater levels and pumping rates and long records of TDS, Mn, and Fe have been used to evaluate the model performance for different runtime stages (training/calibration, testing, and validation). Many of the previous studies have utilized a similar approach [2, 50, 54].The Nash-Sutcliffe Model Efficiency (NSE) is given as(3)NSE=1−∑i=1nGWLio−GWLip2∑i=1nGWLio−GWLavg2,where “GWLio and GWLip)2” denote the observed and predicted groundwater levels, respectively, for the ith data value and “n” is the total number of observed/predicted data points. The NSE values ranging from 0.75 to 1.0 may be termed as “very good,” 0.65 to 0.75 may be taken as “good,” and 0.5 to 0.65 as “satisfactory,” and the values ranging between 0.4 and 0.5 may be expressed as “acceptable” regarding the model’s performance [55, 56].The Mean Square Error (MSE) is given as(4)MSE=1n∑i=1nGWLio−GWLip2.In the case of quality modelling,GWL will be replaced by the values of TDS or the concentrations of Mn and Fe in (3) and (4).


## 3. Results and Discussion

### 3.1. Quantity Modelling

#### 3.1.1. Results of Hydraulic and AI Models

The results of calibration, adjustment of general head boundary condition, and validation are shown in Figures6(a) and 6(b). Figure 6(a) indicates that a value of NSE equal to 0.9 has been achieved. It is quite a reasonable value. It falls in the range of “very good” as per the conditions stated by previous researchers [55, 56]. The simulated versus recorded graph shown in Figure 6(b) is also close to the 1 : 1 line (the 1 : 1 line represents ideal fitting between the simulated and recorded groundwater levels). The NSE equal to 0.87 has been obtained in validation which advocates that the validation results are also acceptable. As mentioned in Section 2.3, one of the highly challenging and original tasks of this research work is the adjustment of the general head boundary condition for the hydraulic model. The results represented in Figures 6(a) and 6(b) show that this task has been executed plausibly. Some of the simulated groundwater depths have been minutely underestimated, whereas others are slightly overestimated. Previous research works by Lyons et al. [57] and Mohanty et al. [58] support these results. Such results are acceptable in groundwater modelling. There are several reasons for this aspect. Firstly, the selection of the hydraulic parameters of the model, including hydraulic conductivity and specific storage, always contributes to the uncertainty in the model results. Secondly, the quality of data in technologically developing countries is not that excellent to obtain very high efficiencies of a groundwater model. Thirdly, the local grid refinement may have a significant impact on parameter estimates, which can ultimately affect the modelling results. Finally, the accuracy of the results also depends upon how sophisticated the optimization technique that has been used in MODFLOW. In short, it is concluded from these results that the traditional methods (physical and numerical models) rely on various input values of parameters and the underlying mechanisms are normally too complicated to grasp. Hence the data-driven approaches should also be tested.Figure 6
MODFLOW results. (a) The values of NSE. (b) Comparison of simulated versus observed groundwater depths and 1 : 1 line (45o line).
(a)(b)The results of ANN are shown in Figures7(a)and 7(b). The hit-and-trial method was used to find the optimal value of neurons in the case of all 5 ANN models (ANN1, ANN2, ANN3, ANN4, and ANN5). The best performance was noticed for ANNs with 10 hidden neurons in both the two-layer architecture and the three layer architecture. The performance for 5 or 15 neurons was not as good as that for 10 neurons. Within an ANN, a neuron is a mathematical function that models the functioning of a biological neuron. There is a lot of literature on the function of neurons. How the behavior of ANN changes with changing the number of neurons is also a topic of reading. However, engineers are interested in obtaining the best solution of an applied problem. Therefore, the exercise performed in this paper is sufficient for solving the groundwater problem at hand. So, the use of 10 neurons has been adopted in the comparison of ANN for the double layer (DL) and the triple layer (TL). This exercise was done for model ANN1 only. The values of NSE equal to 0.97, 0.81, and 0.92 for training, testing, and validation, respectively, in case of TL have been obtained in comparison to the NSE values of 0.95, 0.94, and 0.96 for training, testing, and validation, respectively, in case of TL. The performance of the model with NSE values ranging from 0.75 to 1.0 may be termed as “very good.” The MSE error is also very small in the range of 0.002 to 0.004. Hence the performances of both models are “very good.” However, the ANN with TL performs comparatively better than the ANN with DL. Deep understanding of the behaviors of ANNs with DL and TL needs a lot of reading of the literature on ANN, but, for engineers working in the field of groundwater hydrology, the results of this part of the paper may be sufficient. Similar results have been reported by Stylianoudaki et al. [50] and Almuhaylan et al. [2]. Hence, for further analysis, the ANNs with 10 neurons and TL have been used.Figure 7
ANN1 results. (a) The values of NSE. (b) The values of MSE.
(a)(b)The performances of five ANN models are shown in Figures8(a) and 8(b). According to Figure 8(a), the performance of ANN1 (the Scaled Conjugate Gradient training function) is observed to be comparatively better than those of the other ANN models (ANN2: Levenberg-Marquardt training function, ANN3: Bayesian regularization training function, ANN4: BFGS Quasi-Newton BP training function, and ANN5: resilient backpropagation training function). The Scaled Conjugate Gradient is an efficient training algorithm having highly quick convergence with a very high degree of accuracy. Most probably, comparatively better performance of ANN1 may be due to the reason that a highly effective optimization routine has been used in the Scaled Conjugate Gradient for obtaining the best possible global minimum [2, 57, 58]. However, Figure 8(a) shows that the NSE values range from 0.65 to 0.97 for almost all the ANN models. Hence, the overall performance of all the five models can be categorized as “good” to “very good.” The value of Mean Square Error (Figure 8(b)) ranging from 0.02 to 0.03 further strengthens these findings. Quite similar predictions have been reported by the previous research works [2, 17, 50, 59, 60]. However, most of these publications lack in one way or the other covering a wide range of aspects of groundwater flow, which makes the present research a comprehensive study and an original contribution.Figure 8
ANN model results. (a) The values of NSE. (b) The values of MSE.
(a)(b)Figures9(a), 9(b), and 9(c) show modelling results regarding the quantity of groundwater by various techniques. Overall values of NSE for training, testing, and validation are shown in these figures. It is observed that the performances of ANN and ANFIS models are comparatively better than that of the hydraulic model (MODFLOW). The values of NSE are noticed to be 0.998 for ANFIS1, 0.997 for ANFIS2, 0.997 for ANN1, 0.995 for ANN2, 0.996 for ANN3, 0.995 for ANN4, and 0.995 for ANN5 to 0.997 for ANN and 0.9 for hydraulic model (Figure 9(a)). The predictions by recent publications also support our simulations [12, 13, 15, 16, 50]. Figures 9(b) and 9(c) further confirm this finding. The simulated versus recorded graph is close to the 1 : 1 line in case of ANN and ANFIS [61–63]. As stated above, the 1 : 1 line represents the best matching between the simulated and recorded groundwater levels. In the case of hydraulic model, the simulated groundwater depths in the beginning have been minutely underestimated, whereas others are slightly overestimated (Figure 9(b)). The reason for this has already been explained in Section 3.1.1.Figure 9
Comparison of the ANN, ANFIS, and hydraulic models. (a) Values of NSE. (b) Predicted versus recorded groundwater depths. (c) Comparison between measured and simulated results by various models with respect to time.
(a)(b)(c)


### 3.2. Quality Modelling

#### 3.2.1. Results of Hydraulic and AI Models

The overall results of ANN, ANFIS, and hydraulic models are compared in Figures10(a)–10(g). In quality modelling, the performances of the double-layer and triple-layer ANN models are found to be nearly similar. There is not much difference between the NSE values for the DL and TL architecture (visual comparison is made in Figures 10(a) and 10(b), and exact values of NSE can be seen in Tables 6 and 7). The MSE values given in Figures 10(c) and 10(d) also confirm this finding.Figure 10
Results of the ANN, ANFIS, and hydraulic models. (a) Values of NSE for DL ANN, (b) values of NSE for TL, (c) values of MSE for DL ANN, (d) values of MSE for TL ANN, (e) predicted versus recorded TDS, (f) predicted versus recorded Fe, and (g) predicted versus recorded Mn.
(a)(b)(c)(d)(e)(f)(g)Table 6
NSE values for DL and TL ANN Models.
ModelTL NSE valuesDL NSE valuesTDSFeMnTDSFeMnANN10.990580.755910.859160.989560.769750.84309ANN20.994730.642380.738680.99490.65940.75012ANN30.913110.63840.729180.657390.816140.94077ANN40.938890.730290.840360.940770.746650.73341ANN50.937360.705630.654950.711820.818110.99058ANFIS10.998150.953080.960030.998150.953080.96003ANFIS20.99870.930510.88660.99870.930510.8866HM0.9890.786650.793410.9890.786650.79341Table 7
MSE values for DL and TL ANN Models.
ModelQuality parametersTL NSE valuesDL NSE valuesTDSFeMnTDSFeMnANN10.0002720.0088110.0027630.0003870.0083240.003474ANN20.0001660.0165620.0052080.0001450.0115910.004689ANN30.0025850.0121320.0049860.0002130.0115750.003553ANN40.0018480.0096360.0031390.0016450.0091230.005025ANN50.0018160.0102320.0060520.0029270.010190.00364ANFIS15.24E−051.87E−031.64E−035.24E−051.87E−031.64E−03ANFIS29.01E−050.002730.0022499.01E−050.002730.002249HM0.0012960.0089110.0029630.0012960.010190.00414Regarding five ANN models, the performance of ANN1 is outstanding. The NSE values are 0.99, 0.756, and 0.859 for TDS, Fe, and Mn, respectively, in the case of TL and 0.989, 0.769, and 0.843 for TDS, Fe, and Mn, respectively, in the case of DL architecture of model. However, the performances of the other four ANN models are also acceptable (Figures10(a)–10(g), Figures 11–11(d), and Table 6 and 7). Only five or six values of NSE out of 30 are less than 0.7.Figure 11
Taylor’s Diagram showing correlation, standard deviation, and Root Mean Square Error (RMSE, represented by blue contours) for (a) TDS, (b) Fe, (c) Mn, and (d) groundwater levels.
(a)(b)(c)(d)According to Figures10–10(d), the performances of ANFIS1 and ANFIS 2 are better than that of ANN. The values of NSE range from 0.886 to 0.998, although the performances of ANN and hydraulic models can also be rated as “good” to “excellent.” The values of NSE for ANN1 in the case of TDS, Fe, and Mn are 0.99, 0.756, and 0.859, respectively. The values of NSE for ANN2 in the case of TDS, Fe, and Mn are 0.995, 0.642, and 0.738, respectively; for ANN3, these values are 0.91, 0.638, and 0.729, respectively; for ANN4, these values are 0.939, 0.73, and 0.84; for ANN5 in the case of TDS, Fe, and Mn, these values are 0.937, 0.71, and 0.655 respectively; for ANFIS1 in the case of TDS, Fe, and Mn, these values are 0.998, 0.953, and; for ANFIS2, these values are 0.96, 0.999, 0.93, and 0.887, respectively. However, Figures 10(e)–10(g) and Taylor’s Diagrams (Figures 11(a)–11(d)) show different behavior. Wide range scatters from 1 : 1 line for predicted values indicate that the performances of the models are not that good, as shown by NSE and MSE. Hence the researchers should not always depend upon the NSE and MSE values only to assess the performance of models. In Taylor’s Diagram, the graphs showing predicted versus measured results and some other statistical parameters should also be investigated. This fact is confirmed by Taylor’s Diagram (Figures 11(a)–11(d)). The discrepancy in model results may be due to several reasons. The code selection, impacts of simplifying assumptions in conceptualization, spatial, temporal resolution, and data accuracy are a few parameters that play an important role in determining the reliability and accuracy of the model predictions.The standard deviation values for different parameters are given in Table8. The maximum concentrations of TDS, Fe, and Mn in the groundwater have been found to be 3553 mg/L, 4.29 mg/L, and 1.49 mg/L, respectively. The corresponding minimum values are 456 mg/L, 0.78 mg/L, and 0.19 mg/L, respectively. The standard deviation and the ranges of minimum and maximum values of TDS and concentrations of Fe and Mn show that there is a large variation in amounts of contaminants. Furthermore, the maximum values of TDS and concentrations of Fe and Mn exceeded the WHO standard limits set for drinking, which are 700 mg/L, 0.3, and 0.1 mg/L, respectively, for TDS, Fe, and Mn. Hence the studied groundwater needs treatment before its use for drinking purposes.Table 8
Standard deviation of data from different sources/models.
ModelTDSFeMnMeasured466.31180.6359190.170473ANN1447.12210.5133610.153118ANN2445.26090.5389020.103719ANN3475.06530.3676940.131233ANN4385.52420.4694380.157027ANN5464.10020.4324640.116558ANFIS1459.89860.5928230.161822ANFIS2459.30050.5797130.152165HM442.94130.4741320.161737


### 3.3. Results of Sensitivity Analysis

It was found from sensitivity analysis that hydraulic conductivity is the most crucial hydraulic parameter as compared to the other parameters. Any small change in hydraulic conductivity results in significant changes in the predicted groundwater levels. Hence it is the key parameter in the groundwater quantity modelling and must be estimated with very high accuracy either by pumping tests or by using a powerful optimization technique in the calibration of the model.The sensitivity analysis of contaminant transport model parameters showed that sensitivity from small to large was as follows: hydraulic conductivity, porosity, and dispersion coefficient. The results of sensitivity analysis are quite realistic. The most important parameter in contaminant transport is the dispersion coefficient.


### 3.4. Future Predictions (8 Scenarios)

The future predictions of groundwater depths have been made by MODFLOW for 8 different scenarios, as shown in Figures12(a) and 12(b). The drawdown values up to the year 2070 with respect to the groundwater levels in 2020 are found to be 52, 70, 82, 90, 105, 118, 135, and 66 for scenarios 1 to 8, respectively (Figure 12(a)). The scenarios have developed with respect to the reference pumping rates in 2020. It is observed that even if the pumping rates are not increased as per increase in the population, the drawdown is very high (70 m). The decrease in pumping rates also results in a drawdown of 52 m. This alarming situation not only threatens the water stress but also warns of a very bad pollution situation (Figure 12(b)). The Saq Aquifer is a confined aquifer. The confining cover has another unconfined aquifer called Qassim Aquifer. It has shallow depths of groundwater which is highly contaminated. Extremely high drawdowns in the Saq Aquifer may cause excessively low pressures in the underlying aquifer, permitting polluted water from overlying Qassim Aquifer to travel downwards. This may cause a severe ecological disaster. Hence adaptation of highly effective sustainable planning, development, and management strategies becomes inevitable.Figure 12
Future prediction results of groundwater depths/aquifer depletion for various pumping rate scenarios.
(a)(b)


## 3.1. Quantity Modelling

### 3.1.1. Results of Hydraulic and AI Models

The results of calibration, adjustment of general head boundary condition, and validation are shown in Figures6(a) and 6(b). Figure 6(a) indicates that a value of NSE equal to 0.9 has been achieved. It is quite a reasonable value. It falls in the range of “very good” as per the conditions stated by previous researchers [55, 56]. The simulated versus recorded graph shown in Figure 6(b) is also close to the 1 : 1 line (the 1 : 1 line represents ideal fitting between the simulated and recorded groundwater levels). The NSE equal to 0.87 has been obtained in validation which advocates that the validation results are also acceptable. As mentioned in Section 2.3, one of the highly challenging and original tasks of this research work is the adjustment of the general head boundary condition for the hydraulic model. The results represented in Figures 6(a) and 6(b) show that this task has been executed plausibly. Some of the simulated groundwater depths have been minutely underestimated, whereas others are slightly overestimated. Previous research works by Lyons et al. [57] and Mohanty et al. [58] support these results. Such results are acceptable in groundwater modelling. There are several reasons for this aspect. Firstly, the selection of the hydraulic parameters of the model, including hydraulic conductivity and specific storage, always contributes to the uncertainty in the model results. Secondly, the quality of data in technologically developing countries is not that excellent to obtain very high efficiencies of a groundwater model. Thirdly, the local grid refinement may have a significant impact on parameter estimates, which can ultimately affect the modelling results. Finally, the accuracy of the results also depends upon how sophisticated the optimization technique that has been used in MODFLOW. In short, it is concluded from these results that the traditional methods (physical and numerical models) rely on various input values of parameters and the underlying mechanisms are normally too complicated to grasp. Hence the data-driven approaches should also be tested.Figure 6
MODFLOW results. (a) The values of NSE. (b) Comparison of simulated versus observed groundwater depths and 1 : 1 line (45o line).
(a)(b)The results of ANN are shown in Figures7(a)and 7(b). The hit-and-trial method was used to find the optimal value of neurons in the case of all 5 ANN models (ANN1, ANN2, ANN3, ANN4, and ANN5). The best performance was noticed for ANNs with 10 hidden neurons in both the two-layer architecture and the three layer architecture. The performance for 5 or 15 neurons was not as good as that for 10 neurons. Within an ANN, a neuron is a mathematical function that models the functioning of a biological neuron. There is a lot of literature on the function of neurons. How the behavior of ANN changes with changing the number of neurons is also a topic of reading. However, engineers are interested in obtaining the best solution of an applied problem. Therefore, the exercise performed in this paper is sufficient for solving the groundwater problem at hand. So, the use of 10 neurons has been adopted in the comparison of ANN for the double layer (DL) and the triple layer (TL). This exercise was done for model ANN1 only. The values of NSE equal to 0.97, 0.81, and 0.92 for training, testing, and validation, respectively, in case of TL have been obtained in comparison to the NSE values of 0.95, 0.94, and 0.96 for training, testing, and validation, respectively, in case of TL. The performance of the model with NSE values ranging from 0.75 to 1.0 may be termed as “very good.” The MSE error is also very small in the range of 0.002 to 0.004. Hence the performances of both models are “very good.” However, the ANN with TL performs comparatively better than the ANN with DL. Deep understanding of the behaviors of ANNs with DL and TL needs a lot of reading of the literature on ANN, but, for engineers working in the field of groundwater hydrology, the results of this part of the paper may be sufficient. Similar results have been reported by Stylianoudaki et al. [50] and Almuhaylan et al. [2]. Hence, for further analysis, the ANNs with 10 neurons and TL have been used.Figure 7
ANN1 results. (a) The values of NSE. (b) The values of MSE.
(a)(b)The performances of five ANN models are shown in Figures8(a) and 8(b). According to Figure 8(a), the performance of ANN1 (the Scaled Conjugate Gradient training function) is observed to be comparatively better than those of the other ANN models (ANN2: Levenberg-Marquardt training function, ANN3: Bayesian regularization training function, ANN4: BFGS Quasi-Newton BP training function, and ANN5: resilient backpropagation training function). The Scaled Conjugate Gradient is an efficient training algorithm having highly quick convergence with a very high degree of accuracy. Most probably, comparatively better performance of ANN1 may be due to the reason that a highly effective optimization routine has been used in the Scaled Conjugate Gradient for obtaining the best possible global minimum [2, 57, 58]. However, Figure 8(a) shows that the NSE values range from 0.65 to 0.97 for almost all the ANN models. Hence, the overall performance of all the five models can be categorized as “good” to “very good.” The value of Mean Square Error (Figure 8(b)) ranging from 0.02 to 0.03 further strengthens these findings. Quite similar predictions have been reported by the previous research works [2, 17, 50, 59, 60]. However, most of these publications lack in one way or the other covering a wide range of aspects of groundwater flow, which makes the present research a comprehensive study and an original contribution.Figure 8
ANN model results. (a) The values of NSE. (b) The values of MSE.
(a)(b)Figures9(a), 9(b), and 9(c) show modelling results regarding the quantity of groundwater by various techniques. Overall values of NSE for training, testing, and validation are shown in these figures. It is observed that the performances of ANN and ANFIS models are comparatively better than that of the hydraulic model (MODFLOW). The values of NSE are noticed to be 0.998 for ANFIS1, 0.997 for ANFIS2, 0.997 for ANN1, 0.995 for ANN2, 0.996 for ANN3, 0.995 for ANN4, and 0.995 for ANN5 to 0.997 for ANN and 0.9 for hydraulic model (Figure 9(a)). The predictions by recent publications also support our simulations [12, 13, 15, 16, 50]. Figures 9(b) and 9(c) further confirm this finding. The simulated versus recorded graph is close to the 1 : 1 line in case of ANN and ANFIS [61–63]. As stated above, the 1 : 1 line represents the best matching between the simulated and recorded groundwater levels. In the case of hydraulic model, the simulated groundwater depths in the beginning have been minutely underestimated, whereas others are slightly overestimated (Figure 9(b)). The reason for this has already been explained in Section 3.1.1.Figure 9
Comparison of the ANN, ANFIS, and hydraulic models. (a) Values of NSE. (b) Predicted versus recorded groundwater depths. (c) Comparison between measured and simulated results by various models with respect to time.
(a)(b)(c)


## 3.1.1. Results of Hydraulic and AI Models

The results of calibration, adjustment of general head boundary condition, and validation are shown in Figures6(a) and 6(b). Figure 6(a) indicates that a value of NSE equal to 0.9 has been achieved. It is quite a reasonable value. It falls in the range of “very good” as per the conditions stated by previous researchers [55, 56]. The simulated versus recorded graph shown in Figure 6(b) is also close to the 1 : 1 line (the 1 : 1 line represents ideal fitting between the simulated and recorded groundwater levels). The NSE equal to 0.87 has been obtained in validation which advocates that the validation results are also acceptable. As mentioned in Section 2.3, one of the highly challenging and original tasks of this research work is the adjustment of the general head boundary condition for the hydraulic model. The results represented in Figures 6(a) and 6(b) show that this task has been executed plausibly. Some of the simulated groundwater depths have been minutely underestimated, whereas others are slightly overestimated. Previous research works by Lyons et al. [57] and Mohanty et al. [58] support these results. Such results are acceptable in groundwater modelling. There are several reasons for this aspect. Firstly, the selection of the hydraulic parameters of the model, including hydraulic conductivity and specific storage, always contributes to the uncertainty in the model results. Secondly, the quality of data in technologically developing countries is not that excellent to obtain very high efficiencies of a groundwater model. Thirdly, the local grid refinement may have a significant impact on parameter estimates, which can ultimately affect the modelling results. Finally, the accuracy of the results also depends upon how sophisticated the optimization technique that has been used in MODFLOW. In short, it is concluded from these results that the traditional methods (physical and numerical models) rely on various input values of parameters and the underlying mechanisms are normally too complicated to grasp. Hence the data-driven approaches should also be tested.Figure 6
MODFLOW results. (a) The values of NSE. (b) Comparison of simulated versus observed groundwater depths and 1 : 1 line (45o line).
(a)(b)The results of ANN are shown in Figures7(a)and 7(b). The hit-and-trial method was used to find the optimal value of neurons in the case of all 5 ANN models (ANN1, ANN2, ANN3, ANN4, and ANN5). The best performance was noticed for ANNs with 10 hidden neurons in both the two-layer architecture and the three layer architecture. The performance for 5 or 15 neurons was not as good as that for 10 neurons. Within an ANN, a neuron is a mathematical function that models the functioning of a biological neuron. There is a lot of literature on the function of neurons. How the behavior of ANN changes with changing the number of neurons is also a topic of reading. However, engineers are interested in obtaining the best solution of an applied problem. Therefore, the exercise performed in this paper is sufficient for solving the groundwater problem at hand. So, the use of 10 neurons has been adopted in the comparison of ANN for the double layer (DL) and the triple layer (TL). This exercise was done for model ANN1 only. The values of NSE equal to 0.97, 0.81, and 0.92 for training, testing, and validation, respectively, in case of TL have been obtained in comparison to the NSE values of 0.95, 0.94, and 0.96 for training, testing, and validation, respectively, in case of TL. The performance of the model with NSE values ranging from 0.75 to 1.0 may be termed as “very good.” The MSE error is also very small in the range of 0.002 to 0.004. Hence the performances of both models are “very good.” However, the ANN with TL performs comparatively better than the ANN with DL. Deep understanding of the behaviors of ANNs with DL and TL needs a lot of reading of the literature on ANN, but, for engineers working in the field of groundwater hydrology, the results of this part of the paper may be sufficient. Similar results have been reported by Stylianoudaki et al. [50] and Almuhaylan et al. [2]. Hence, for further analysis, the ANNs with 10 neurons and TL have been used.Figure 7
ANN1 results. (a) The values of NSE. (b) The values of MSE.
(a)(b)The performances of five ANN models are shown in Figures8(a) and 8(b). According to Figure 8(a), the performance of ANN1 (the Scaled Conjugate Gradient training function) is observed to be comparatively better than those of the other ANN models (ANN2: Levenberg-Marquardt training function, ANN3: Bayesian regularization training function, ANN4: BFGS Quasi-Newton BP training function, and ANN5: resilient backpropagation training function). The Scaled Conjugate Gradient is an efficient training algorithm having highly quick convergence with a very high degree of accuracy. Most probably, comparatively better performance of ANN1 may be due to the reason that a highly effective optimization routine has been used in the Scaled Conjugate Gradient for obtaining the best possible global minimum [2, 57, 58]. However, Figure 8(a) shows that the NSE values range from 0.65 to 0.97 for almost all the ANN models. Hence, the overall performance of all the five models can be categorized as “good” to “very good.” The value of Mean Square Error (Figure 8(b)) ranging from 0.02 to 0.03 further strengthens these findings. Quite similar predictions have been reported by the previous research works [2, 17, 50, 59, 60]. However, most of these publications lack in one way or the other covering a wide range of aspects of groundwater flow, which makes the present research a comprehensive study and an original contribution.Figure 8
ANN model results. (a) The values of NSE. (b) The values of MSE.
(a)(b)Figures9(a), 9(b), and 9(c) show modelling results regarding the quantity of groundwater by various techniques. Overall values of NSE for training, testing, and validation are shown in these figures. It is observed that the performances of ANN and ANFIS models are comparatively better than that of the hydraulic model (MODFLOW). The values of NSE are noticed to be 0.998 for ANFIS1, 0.997 for ANFIS2, 0.997 for ANN1, 0.995 for ANN2, 0.996 for ANN3, 0.995 for ANN4, and 0.995 for ANN5 to 0.997 for ANN and 0.9 for hydraulic model (Figure 9(a)). The predictions by recent publications also support our simulations [12, 13, 15, 16, 50]. Figures 9(b) and 9(c) further confirm this finding. The simulated versus recorded graph is close to the 1 : 1 line in case of ANN and ANFIS [61–63]. As stated above, the 1 : 1 line represents the best matching between the simulated and recorded groundwater levels. In the case of hydraulic model, the simulated groundwater depths in the beginning have been minutely underestimated, whereas others are slightly overestimated (Figure 9(b)). The reason for this has already been explained in Section 3.1.1.Figure 9
Comparison of the ANN, ANFIS, and hydraulic models. (a) Values of NSE. (b) Predicted versus recorded groundwater depths. (c) Comparison between measured and simulated results by various models with respect to time.
(a)(b)(c)


## 3.2. Quality Modelling

### 3.2.1. Results of Hydraulic and AI Models

The overall results of ANN, ANFIS, and hydraulic models are compared in Figures10(a)–10(g). In quality modelling, the performances of the double-layer and triple-layer ANN models are found to be nearly similar. There is not much difference between the NSE values for the DL and TL architecture (visual comparison is made in Figures 10(a) and 10(b), and exact values of NSE can be seen in Tables 6 and 7). The MSE values given in Figures 10(c) and 10(d) also confirm this finding.Figure 10
Results of the ANN, ANFIS, and hydraulic models. (a) Values of NSE for DL ANN, (b) values of NSE for TL, (c) values of MSE for DL ANN, (d) values of MSE for TL ANN, (e) predicted versus recorded TDS, (f) predicted versus recorded Fe, and (g) predicted versus recorded Mn.
(a)(b)(c)(d)(e)(f)(g)Table 6
NSE values for DL and TL ANN Models.
ModelTL NSE valuesDL NSE valuesTDSFeMnTDSFeMnANN10.990580.755910.859160.989560.769750.84309ANN20.994730.642380.738680.99490.65940.75012ANN30.913110.63840.729180.657390.816140.94077ANN40.938890.730290.840360.940770.746650.73341ANN50.937360.705630.654950.711820.818110.99058ANFIS10.998150.953080.960030.998150.953080.96003ANFIS20.99870.930510.88660.99870.930510.8866HM0.9890.786650.793410.9890.786650.79341Table 7
MSE values for DL and TL ANN Models.
ModelQuality parametersTL NSE valuesDL NSE valuesTDSFeMnTDSFeMnANN10.0002720.0088110.0027630.0003870.0083240.003474ANN20.0001660.0165620.0052080.0001450.0115910.004689ANN30.0025850.0121320.0049860.0002130.0115750.003553ANN40.0018480.0096360.0031390.0016450.0091230.005025ANN50.0018160.0102320.0060520.0029270.010190.00364ANFIS15.24E−051.87E−031.64E−035.24E−051.87E−031.64E−03ANFIS29.01E−050.002730.0022499.01E−050.002730.002249HM0.0012960.0089110.0029630.0012960.010190.00414Regarding five ANN models, the performance of ANN1 is outstanding. The NSE values are 0.99, 0.756, and 0.859 for TDS, Fe, and Mn, respectively, in the case of TL and 0.989, 0.769, and 0.843 for TDS, Fe, and Mn, respectively, in the case of DL architecture of model. However, the performances of the other four ANN models are also acceptable (Figures10(a)–10(g), Figures 11–11(d), and Table 6 and 7). Only five or six values of NSE out of 30 are less than 0.7.Figure 11
Taylor’s Diagram showing correlation, standard deviation, and Root Mean Square Error (RMSE, represented by blue contours) for (a) TDS, (b) Fe, (c) Mn, and (d) groundwater levels.
(a)(b)(c)(d)According to Figures10–10(d), the performances of ANFIS1 and ANFIS 2 are better than that of ANN. The values of NSE range from 0.886 to 0.998, although the performances of ANN and hydraulic models can also be rated as “good” to “excellent.” The values of NSE for ANN1 in the case of TDS, Fe, and Mn are 0.99, 0.756, and 0.859, respectively. The values of NSE for ANN2 in the case of TDS, Fe, and Mn are 0.995, 0.642, and 0.738, respectively; for ANN3, these values are 0.91, 0.638, and 0.729, respectively; for ANN4, these values are 0.939, 0.73, and 0.84; for ANN5 in the case of TDS, Fe, and Mn, these values are 0.937, 0.71, and 0.655 respectively; for ANFIS1 in the case of TDS, Fe, and Mn, these values are 0.998, 0.953, and; for ANFIS2, these values are 0.96, 0.999, 0.93, and 0.887, respectively. However, Figures 10(e)–10(g) and Taylor’s Diagrams (Figures 11(a)–11(d)) show different behavior. Wide range scatters from 1 : 1 line for predicted values indicate that the performances of the models are not that good, as shown by NSE and MSE. Hence the researchers should not always depend upon the NSE and MSE values only to assess the performance of models. In Taylor’s Diagram, the graphs showing predicted versus measured results and some other statistical parameters should also be investigated. This fact is confirmed by Taylor’s Diagram (Figures 11(a)–11(d)). The discrepancy in model results may be due to several reasons. The code selection, impacts of simplifying assumptions in conceptualization, spatial, temporal resolution, and data accuracy are a few parameters that play an important role in determining the reliability and accuracy of the model predictions.The standard deviation values for different parameters are given in Table8. The maximum concentrations of TDS, Fe, and Mn in the groundwater have been found to be 3553 mg/L, 4.29 mg/L, and 1.49 mg/L, respectively. The corresponding minimum values are 456 mg/L, 0.78 mg/L, and 0.19 mg/L, respectively. The standard deviation and the ranges of minimum and maximum values of TDS and concentrations of Fe and Mn show that there is a large variation in amounts of contaminants. Furthermore, the maximum values of TDS and concentrations of Fe and Mn exceeded the WHO standard limits set for drinking, which are 700 mg/L, 0.3, and 0.1 mg/L, respectively, for TDS, Fe, and Mn. Hence the studied groundwater needs treatment before its use for drinking purposes.Table 8
Standard deviation of data from different sources/models.
ModelTDSFeMnMeasured466.31180.6359190.170473ANN1447.12210.5133610.153118ANN2445.26090.5389020.103719ANN3475.06530.3676940.131233ANN4385.52420.4694380.157027ANN5464.10020.4324640.116558ANFIS1459.89860.5928230.161822ANFIS2459.30050.5797130.152165HM442.94130.4741320.161737


## 3.2.1. Results of Hydraulic and AI Models

The overall results of ANN, ANFIS, and hydraulic models are compared in Figures10(a)–10(g). In quality modelling, the performances of the double-layer and triple-layer ANN models are found to be nearly similar. There is not much difference between the NSE values for the DL and TL architecture (visual comparison is made in Figures 10(a) and 10(b), and exact values of NSE can be seen in Tables 6 and 7). The MSE values given in Figures 10(c) and 10(d) also confirm this finding.Figure 10
Results of the ANN, ANFIS, and hydraulic models. (a) Values of NSE for DL ANN, (b) values of NSE for TL, (c) values of MSE for DL ANN, (d) values of MSE for TL ANN, (e) predicted versus recorded TDS, (f) predicted versus recorded Fe, and (g) predicted versus recorded Mn.
(a)(b)(c)(d)(e)(f)(g)Table 6
NSE values for DL and TL ANN Models.
ModelTL NSE valuesDL NSE valuesTDSFeMnTDSFeMnANN10.990580.755910.859160.989560.769750.84309ANN20.994730.642380.738680.99490.65940.75012ANN30.913110.63840.729180.657390.816140.94077ANN40.938890.730290.840360.940770.746650.73341ANN50.937360.705630.654950.711820.818110.99058ANFIS10.998150.953080.960030.998150.953080.96003ANFIS20.99870.930510.88660.99870.930510.8866HM0.9890.786650.793410.9890.786650.79341Table 7
MSE values for DL and TL ANN Models.
ModelQuality parametersTL NSE valuesDL NSE valuesTDSFeMnTDSFeMnANN10.0002720.0088110.0027630.0003870.0083240.003474ANN20.0001660.0165620.0052080.0001450.0115910.004689ANN30.0025850.0121320.0049860.0002130.0115750.003553ANN40.0018480.0096360.0031390.0016450.0091230.005025ANN50.0018160.0102320.0060520.0029270.010190.00364ANFIS15.24E−051.87E−031.64E−035.24E−051.87E−031.64E−03ANFIS29.01E−050.002730.0022499.01E−050.002730.002249HM0.0012960.0089110.0029630.0012960.010190.00414Regarding five ANN models, the performance of ANN1 is outstanding. The NSE values are 0.99, 0.756, and 0.859 for TDS, Fe, and Mn, respectively, in the case of TL and 0.989, 0.769, and 0.843 for TDS, Fe, and Mn, respectively, in the case of DL architecture of model. However, the performances of the other four ANN models are also acceptable (Figures10(a)–10(g), Figures 11–11(d), and Table 6 and 7). Only five or six values of NSE out of 30 are less than 0.7.Figure 11
Taylor’s Diagram showing correlation, standard deviation, and Root Mean Square Error (RMSE, represented by blue contours) for (a) TDS, (b) Fe, (c) Mn, and (d) groundwater levels.
(a)(b)(c)(d)According to Figures10–10(d), the performances of ANFIS1 and ANFIS 2 are better than that of ANN. The values of NSE range from 0.886 to 0.998, although the performances of ANN and hydraulic models can also be rated as “good” to “excellent.” The values of NSE for ANN1 in the case of TDS, Fe, and Mn are 0.99, 0.756, and 0.859, respectively. The values of NSE for ANN2 in the case of TDS, Fe, and Mn are 0.995, 0.642, and 0.738, respectively; for ANN3, these values are 0.91, 0.638, and 0.729, respectively; for ANN4, these values are 0.939, 0.73, and 0.84; for ANN5 in the case of TDS, Fe, and Mn, these values are 0.937, 0.71, and 0.655 respectively; for ANFIS1 in the case of TDS, Fe, and Mn, these values are 0.998, 0.953, and; for ANFIS2, these values are 0.96, 0.999, 0.93, and 0.887, respectively. However, Figures 10(e)–10(g) and Taylor’s Diagrams (Figures 11(a)–11(d)) show different behavior. Wide range scatters from 1 : 1 line for predicted values indicate that the performances of the models are not that good, as shown by NSE and MSE. Hence the researchers should not always depend upon the NSE and MSE values only to assess the performance of models. In Taylor’s Diagram, the graphs showing predicted versus measured results and some other statistical parameters should also be investigated. This fact is confirmed by Taylor’s Diagram (Figures 11(a)–11(d)). The discrepancy in model results may be due to several reasons. The code selection, impacts of simplifying assumptions in conceptualization, spatial, temporal resolution, and data accuracy are a few parameters that play an important role in determining the reliability and accuracy of the model predictions.The standard deviation values for different parameters are given in Table8. The maximum concentrations of TDS, Fe, and Mn in the groundwater have been found to be 3553 mg/L, 4.29 mg/L, and 1.49 mg/L, respectively. The corresponding minimum values are 456 mg/L, 0.78 mg/L, and 0.19 mg/L, respectively. The standard deviation and the ranges of minimum and maximum values of TDS and concentrations of Fe and Mn show that there is a large variation in amounts of contaminants. Furthermore, the maximum values of TDS and concentrations of Fe and Mn exceeded the WHO standard limits set for drinking, which are 700 mg/L, 0.3, and 0.1 mg/L, respectively, for TDS, Fe, and Mn. Hence the studied groundwater needs treatment before its use for drinking purposes.Table 8
Standard deviation of data from different sources/models.
ModelTDSFeMnMeasured466.31180.6359190.170473ANN1447.12210.5133610.153118ANN2445.26090.5389020.103719ANN3475.06530.3676940.131233ANN4385.52420.4694380.157027ANN5464.10020.4324640.116558ANFIS1459.89860.5928230.161822ANFIS2459.30050.5797130.152165HM442.94130.4741320.161737


## 3.3. Results of Sensitivity Analysis

It was found from sensitivity analysis that hydraulic conductivity is the most crucial hydraulic parameter as compared to the other parameters. Any small change in hydraulic conductivity results in significant changes in the predicted groundwater levels. Hence it is the key parameter in the groundwater quantity modelling and must be estimated with very high accuracy either by pumping tests or by using a powerful optimization technique in the calibration of the model.The sensitivity analysis of contaminant transport model parameters showed that sensitivity from small to large was as follows: hydraulic conductivity, porosity, and dispersion coefficient. The results of sensitivity analysis are quite realistic. The most important parameter in contaminant transport is the dispersion coefficient.


## 3.4. Future Predictions (8 Scenarios)

The future predictions of groundwater depths have been made by MODFLOW for 8 different scenarios, as shown in Figures12(a) and 12(b). The drawdown values up to the year 2070 with respect to the groundwater levels in 2020 are found to be 52, 70, 82, 90, 105, 118, 135, and 66 for scenarios 1 to 8, respectively (Figure 12(a)). The scenarios have developed with respect to the reference pumping rates in 2020. It is observed that even if the pumping rates are not increased as per increase in the population, the drawdown is very high (70 m). The decrease in pumping rates also results in a drawdown of 52 m. This alarming situation not only threatens the water stress but also warns of a very bad pollution situation (Figure 12(b)). The Saq Aquifer is a confined aquifer. The confining cover has another unconfined aquifer called Qassim Aquifer. It has shallow depths of groundwater which is highly contaminated. Extremely high drawdowns in the Saq Aquifer may cause excessively low pressures in the underlying aquifer, permitting polluted water from overlying Qassim Aquifer to travel downwards. This may cause a severe ecological disaster. Hence adaptation of highly effective sustainable planning, development, and management strategies becomes inevitable.Figure 12
Future prediction results of groundwater depths/aquifer depletion for various pumping rate scenarios.
(a)(b)


## 4. Conclusions and Recommendations

The MODFLOW, a well-known hydraulic model, has been applied to forecast the variations in the groundwater quantity and quality of the Saq Aquifer in the Qassim Region. AI models have also efficiently simulated the quantity and quality of groundwater in this aquifer. Different architectures of ANN with triple and double layers having 5, 10, and 15 hidden neurons in every layer have been investigated. Five different categories of training functions are examined in ANN models. Furthermore, two versions of ANFIS have been utilized successfully to predict the depletion of aquifer with respect to both quantity and quality. Three main quality parameters, TDS, Fe, and Mn, have been investigated.The study concluded that both versions of ANFIS (ANFIS1 and ANFIS2) can efficiently predict groundwater levels and contaminants (TDS, Fe, and Mn) with NSE up to 0.99. The ANN model provides the best results with 10 hidden neurons in each hidden layer. The performance of ANN with the architecture having three layers and 10 neurons in each layer is better than the one with two layers in case of the quality and quantity modelling of groundwater. The Scaled Conjugate Gradient training function in ANN has comparatively better performance compared to the Levenberg-Marquardt, Bayesian regularization, BFGS Quasi-Newton BP, and resilient backpropagation training functions for predicting groundwater levels and contaminants. The performance of the hydraulic model is good, but, for the given set of data of pumping wells and contaminants, its performance is not as good as that of the ANN or the ANFIS models. However, it is the most robust and reliable model based on the laws of physics. The hydraulic modelling is comparatively more demanding than the ANN and ANFIS models for both the contaminant and flow transport.There will be extremely excessive drawdowns in groundwater levels of the Saq Aquifer (ranging from 70 to 135 m with respect to the reference water levels in 2020) for various scenarios of pumping rates. Even with the constant rate of pumping without looking at the needs of the increasing population in the region, the drawdown is 70 m for the coming 50 years. The concentrations of contaminants (TDS, Fe, and Mn) are also increasing significantly. There may be a major ecological disaster if preventative actions are not adopted. However, the constant rate of pumping may result in a comparatively increased life of the Saq Aquifer as compared to the other scenarios of pumping. The drawdown will reduce by 33% over the next 50 years’ period compared to the pumping at increased rates to meet the water needs of a growing population. There is a negligible amount of recharge to the Saq Aquifer, which may fade away, causing severe scarcity of water in the future. Thereby, some alternative sources of water should be established to fulfil the goals of the Kingdom’s Vision 2030.The results and models developed in this paper may be very useful in obtaining the pumping rates for an environment-friendly future. For instance, by keeping pumping rates constant (no increase in pumping), there may be a substantial increase in the Saq Aquifer’s life. A 33% decrease in drawdown may be accomplished in the coming 50 years by implementing groundwater preservation strategies compared to the situation of uninterrupted increase in pumping rates which otherwise may be necessary to meet the water requirements of the promptly rising population. Results of this research will assist in the planning, development, and management of stressed water resources in the Gulf Region, Saudi Arabia, and the arid regions with analogous water scarcity situations.An extensive research work to preserve the existing sources and develop new sources of water is recommended. Highly efficient technology, extremely accurate and sufficient data, and demarcation of the location of outcrops are of utmost importance for the study area.



---

*Source: 5910989-2022-05-16.xml*