Application of Probabilistic and Machine Learning Models for Groundwater Potentiality Mapping in Damghan Sedimentary Plain, Iran

: Groundwater is one of the most important natural resources, as it regulates the earth’s hydrological system. The Damghan sedimentary plain area, located in the region of a semi-arid climate of Iran, has very critical conditions of groundwater due to massive pressure on it and is in need of robust models for identifying the groundwater potential zones (GWPZ). The main goal of the current research is to prepare a groundwater potentiality map (GWPM) considering the probabilistic, machine learning, data mining, and multi-criteria decision

. Location of the study area in Iran and Semnan province and location of training and validations wells in the study area.

Methodology
For assessing the groundwater potentiality (GWP), some spatial and non-spatial data have been gathered to prepare different datasets for modeling and validation of results. The data consists of two, i.e., primary and secondary, data. Primary data are pumping tests and yield measurements in the field. The secondary data are topographical map (scale 1:50000), lithological map (scale 1:100000), Sentinel 2A, Phased Array type L-band synthetic aperture radar (PALSAR) digital elevation model (DEM), rainfall of different metrological station of last 30 years, well location data from Water Resource Management, Iran, soil data from soil department of Iran. Thematic maps of all the data were extracted and analyzed by the RS and GIS. The present work methodologically consists of four phases (Figure 2) including; (1) preparation of groundwater inventory database thematic data layers of the groundwater conditioning factors including elevation, slope, aspect, CI, rainfall, lithology, soil type, LU/LC, Dd, distance to river, distance to fault, distance to road, NDVI, TPI), TWI, and SPI; (2) multicollinearity assessment of the effective groundwater determining factors (GWDFs); (3) application of models and preparation of GWPMs. The GWPMs were classified according to the four The sedimentary plain of Damghan is located in arid and semi-arid regions and facing the problem of water supply such as other arid regions. The main source of freshwater is sub-surface water storage and undergoes the problems of over pumping and lowering of groundwater. In the recent decade, due to excessive groundwater exploitation for irrigation and industrial purposes combining with the decreasing amount of rainfall, the water table is coming down rapidly. Therefore, immediate planning is needed to conserve the groundwater [30]. In this respect, the delineation of GWPZ is essential for proper planning and sustainable management of water.

Methodology
For assessing the groundwater potentiality (GWP), some spatial and non-spatial data have been gathered to prepare different datasets for modeling and validation of results. The data consists of two, i.e., primary and secondary, data. Primary data are pumping tests and yield measurements in the field. The secondary data are topographical map (scale 1:50,000), lithological map (scale Remote Sens. 2019, 11, 3015 5 of 35 1:100,000), Sentinel 2A, Phased Array type L-band synthetic aperture radar (PALSAR) digital elevation model (DEM), rainfall of different metrological station of last 30 years, well location data from Water Resource Management, Iran, soil data from soil department of Iran. Thematic maps of all the data were extracted and analyzed by the RS and GIS. The present work methodologically consists of four phases (Figure 2) including; (1) preparation of groundwater inventory database thematic data layers of the groundwater conditioning factors including elevation, slope, aspect, CI, rainfall, lithology, soil type, LU/LC, Dd, distance to river, distance to fault, distance to road, NDVI, TPI), TWI, and SPI; (2) multicollinearity assessment of the effective groundwater determining factors (GWDFs); (3) application of models and preparation of GWPMs. The GWPMs were classified according to the four classification methods, namely quantile, natural breaks, equal interval, and geometrical interval, into four different groundwater susceptibility classes, including low, medium, high, and very high. By comparing the results of each classification method and the distribution of training and validation wells on the high and very high groundwater susceptibility classes, it was found that the natural break classification method gave the most accurate distribution. This agrees with the findings by Arabameri et al. [32], in that natural break method is a good classifier in susceptibility mapping; and (4), evaluation of the models performances using area under receiver operating characteristics (AUROC) curve, sensitivity (SE), specificity (SP), accuracy (AC), mean absolute error (MAE), root mean square error (RMSE) and seed cell area index (SCAI) methods.

Groundwater Inventory Map (GWIM)
The groundwater inventory database is of a key role in groundwater potentiality mapping. An inventory map is a target variable for any spatial modeling [32]. The well inventory database was prepared after extensive field visit with a hand GPS (global positioning system), and yield data were collected from the Department of Water Resources Management, Iran. Groundwater wells, with high yield of ≥11 m 3 h−1 by pumping test analysis, have been considered for the GWPM. As a result, 80 wells have recognized in the study area. 56 wells (70%) of this dataset, were randomly selected to produce the GWPM models [32], whereas the remaining 24 (30%) wells were considered for validation of GWPMs [11]. The training and testing wells locations have been mentioned in Figure 1.

Groundwater Determining Factors (GWDFs)
The different geo-environmental components play a crucial role in determining the status of groundwater. GWPM represents the association between GWDFs and well locations [21,22]. For the GWPM, 16 GWDFs have been selected including elevation, slope, aspect, CI, rainfall, lithology, soil type, LU/LC, NDVI, Dd, distance to the river, distance to fault, distance to road, TWI, SPI and TPI Remote Sens. 2019, 11, 3015 7 of 35 ( Figure 3a-p). The PALSAR DEM (12.5 m resolution) downloaded from the Alaska Satellite Facility (ASF) Distributed Active Archive Center (DAAC). In this study, PALSAR DEM was used to extract the topographical, hydrological factors such as elevation, slope, aspect, CI, drainage, TWI, SPI, and TPI. The slope, aspect, and elevation are the major topographic components, used to determine the groundwater potentiality, erosion probability, etc. [21]. The DEM has been used as the elevation dataset ( Figure 3a). The altitudinal fluctuation controls climatic conditions and helps to induce various vegetation types and soil development [33]. The slope data layer has been derived from PALSAR DEM by spatial analysis in the GIS environment (Figure 3b). In the same way, the aspect map has also been extracted from PALSAR DEM imagery using a spatial analysis tool (Figure 3d). CI is an important terrain factor that demonstrates the arrangement of relief as a set of channels and ridges. It is developed by Kiss [34]. The convergence index (CI) has been calculated using Equation (1).
where θ indicates the average angle between the aspect of adjacent cells and the direction to the central cell. The CI value ranges from −100 to +100 (Figure 3c). The rainfall map was prepared by the kriging method considering the last 10-year annual rainfall of different stations ( Figure 3e). The drainage was extracted from the topographical map and PALSAR DEM imagery. The Dd was computed based on Horton's morphometric formula (Equation (2)).
where Lu means the total length of all orders streams, A is the area in square kilometer. Finally, the spatial data layer of the Dd has been built using the IDW interpolation method in the GIS environment ( Figure 3f). The fault layer has been taken out from Landsat 7 imagery in ENVI software. The road network has been taken off from the topographical map and Google Earth imagery. The distance to river, fault, road data layers have been built using the Euclidian distance buffering tool and expressed in km (Figure 3g,h,i) [11]. The lithological information for the study area was gathered from the geological department of Iran [29]. The lithology map has been prepared by the digitization process ( Figure 3j). Geologically, the region is composed of nine geological segments, namely A, B, C, D, E, F, G, H, and I ( Figure 3j and Table 1). Soil data was collected from the soil department of Iranian and with the help of the digitized process, the thematic dataset of soil has been produced ( Figure 3k). The LU/LC map has been produced from Sentinel 2A satellite image (12/08/2017) of 10 m, 20 m, and 60 m spatial resolution for each band using the supervised image classification method ( Figure 3l). NDVI has been computed from satellite image (Figure 3m) using Equation (3): where NIR is the near-infrared band or band 8 and red band or band 4. The TWI directly affects the topographic conditions, which control the hydrological process. TWI is the function of slope and the upstream area per unit width orthogonal to the direction of flow [35]. TWI plays a major role in the spatial heterogeneity of hydrological conditions such as soil moisture, underwater flow and slope steady-state [32]. The TWI has been introduced by Beven and Kirkby [36]. TWI is calculated from Equation (4): where AS represents the cumulative area of the catchment (m 2 m −1 ) and β is the slope gradient (degrees). The TWI value ranges from 1.11 to 21.54 ( Figure 3n). The SPI is a calculation of water flow erosive power based on the assumption that discharge is commensurate with a given catchment area [37]. One of the most important factors in controlling slope erosion processes is SPI. Regions with high stream power have high erosion potentiality [38]. From Equation (5), SPI has been calculated: where A S is the upstream contributing area and β is slope gradient (in degrees). The spatial allocation of SPI ranges from 6.27 to 24.44 ( Figure 3o) in the research area. TPI is defined as the difference between the middle point elevation (Z 0 ) and the average elevation (Z) in a predetermined radius around it (R) [39]: The TPI has positive and negative value; a positive value demonstrates that the midpoint is located at a higher place than its average while a negative value indicates a lower place than the average. The TPI range depends not only on variations in altitude but also on landscape units (R) [40]. Where large R values mainly depend on the main units of landscape, and small R values show up lower valleys such as small valleys and ridges. The TPI value ranges from 12.16 to 14.67 in this plain ( Figure 3p). Spatial resolutions of the selected GWDFs are not the same. For preparing the groundwater potential maps of the study area the resolution of PALSAR DEM, i.e., 12.5 m* 12.5 m has been selected as the base scale and all the GWCFs of which scale are greater or lesser than the PALSAR DEM have been resembled into a 12.  Table 2.

Weight of Evidence (WoE) Model
The WoE model is the main Bayesian probability system model in linear logic form and uses non-conditional and conditional probabilities [41]. WoE reveals the spatial association between dependent variable, i.e., well locations and independent variables, i.e., GWDFs. The weight of each class has been assigned by this method using the following equations (Equations (8)- (14)) [42]: where P(B|D) is the conditional probability of B occurring given the presence of D, B is the datasets of GWDFs related to the presence of groundwater well and B indicates the groundwater is absent of the datasets of groundwater conditioning factors. D indicates the presence of well while D stand for the absence of a well, and P is the probability. Whereas, W + i is a positive weight of GWDFs for groundwater occurrence. Conversely, W − i is a negative weight with respect to the absence of groundwater well (unfavorable factors). WoE computation has been started by the pixels counting process between groundwater well locations and GWDFs. The weighted GWDFs factors have been summed up in the raster calculation to generate the single layer of GWPM in the GIS environment using the following Equation (15 RF is the non-parametric multivariate model [43], which can be used for the analysis of regression and classification and variable selections. RF model creates thousands of trees, forming a 'forest' based on the decision rule. Each tree in the RF model depends on a sample of bootstrapped of data using a CART process with a random subset of variables selected at each node. The final decision of the class membership and model (output) has determined according to the majority priority of all decision trees [44]. The trees' ensembles would have performed much better than a single tree. It is important to know that the program can be run by a large number of trees with taking large and too many computational requirements [45]. RF is a very reliable and flexible ensemble classifier, which depends upon the decision trees, that have so many attractive performances such a minimum costly, minimum tendencies for overfitting and also capability of the work with very high dimensional data [46]. The RF model is also a very fast machine learning solution, allowing a highly accurate classification with internal unbiased generalizability estimation during the process of forest construction [47]. The basic merits of RF model arise when the program proceeds including (i) no need of any assumptions regarding the data distribution, (ii) no overfitting problem, (iii) in case of single tree, a low correlation estimated, while the diversity of forest increases the usages of a number of factors, (iv) helps to estimate negative or error using 'out-of-bag'(OOB) data, (v) averages a large number of trees, resulting the low bias and low variance, subsequently, (f), resulting in the excellent prediction for performances [43,48]. Besides, the numerical and categorical data can be incorporated in the RF model. Using the OOB error-index, the variance and covariance between the grids cells can be estimated [48]. The predication values of this model are estimated by the huge amount of decision trees [43]. The presence and absence groundwater wells among the GWDFs can easily be estimated by RF model. In this algorithm, the mean decrease accuracy and Gini are estimated by the RF model to analysis the variable importance of the GWDFs [30]. This algorithm calculates untouched the proper count, the amount of correct classification applying the data out-of-bag as its test sample. In the out-of-bag instances, the values of the attributes are then randomly permuted. A new set of data will then be checked for proper classification. The average of this number is the raw importance score for the specific attribute over all trees in the forest. Therefore, factors importance in RF model is computed for variable Yi by out of bag error (OOB) [49]. Factors importance of Yi can be calculated using the Equation (16): where ntree stands for the number of trees, VImp(Yi) denotes variable importance for variable Yj, errOOB t is an error when all the factors are included, and errOOB j t denotes an error after the removal of the variable j. The Gini index was used to measure the variable significance based on the number of times that variable is picked by all trees [47,50]. In this study, the 'randomForest' package in R program has been installed and run the RF model for estimating the GWPZ [51]. Finally, the RF model-based GWPM has been produced in the GIS environment.

Binary Logistic Regression (BLR)
BLR model is the most common statistical model which considers both dichotomous and continuous variables. However, the practical dependent variable must be in binary form, i.e., 1 and 0. Where, 0 represents the absence of groundwater well and 1 stands for the presence of the groundwater well [37,52]. For GWPM, it corresponds with the Bernoulli method, which determines the high groundwater potentiality over space depending on the Bernoulli probability [32]. The main target of the BLR analysis is to chalk out the correct and appropriate prediction of samples and probe the correlation between a dependent variable with a set of independent variables [32]. Among the different methods of the regressions, the BLR is fitting a logistic curve or function concerning data. As a result, BLR estimates the values that vary from 0 to 1, while 1 means the presence of groundwater well, conversely 0 means the chances of occurrence of groundwater well is nil. In this method, the target value is calculated using Equation (17): where Logit is the link function, P is the probability of occurrence of groundwater well (y), p = 1 − p are the odds of groundwater occurrence (or probability of presence divided by the probability of absence) the, C0 is the model intercept and (C1, . . . , Cn) are the regression coefficients for each GWCF (X1, . . . , Xn) [32]. In this contribution, the BLR model has been applied in R by using the 'glm' function based on the "stats" package [32]. In this study, the random point's values have been extracted from each variable of GWCFs for presence and absence condition of the groundwater. Finally, GWPM by BLR model has been produced in GIS with regarding the prediction database.

Technique for Order Preference by Similarity to Ideal Solution (TOPSIS)
The TOPSIS method was introduced by Hwang and Yoon [53]. Presently, it is an important multi-criteria decision approach among the several MCDA processes utilized for water management practices [54]. This approach relies on the premise that the best alternative solution should have the shortest distance Euclidean from the positive ideal solution and the longest distance from the negative solution [32]. The ideal solution in GWPM is the best model to distinguish between groundwater well presences and absences. TOPSIS is a method that cannot understand categorical properties [55]. That is why the AHP has been used to assign the weight to each GWDF. The computation of the entire TOPSIS model was carried out step by step as follows: Step 1: Preparation of a decision matrix with m criteria and n alternatives using Equation (18): Step 2: Normalization of decision matrix using Equation (19): where i = 1, . . . , m; j = 1 . . . n.
Step 3: Determine the weight of criteria using the AHP model: AHP, first implemented by Saaty [56] is one of the most comprehensive MDCA approaches. This method assists decision-makers in receiving quantitative and qualitative parameters. The pairwise comparisons help for the judgment and computation of GWCFs [57]. After preparing paired comparisons, the resulting paired comparison matrix is normalized using Equation (20) and then the final weight (W i ) of each parameter is obtained using Equation (21).
One of the advantages of this method is to show the inconsistency [58]. To assess the degree of weighting precision, an Inconsistency Index is used. The consistency test shows how much trust can be put in the priorities of a matrix. If this value is >0.1, this means that it is not consistent with the specified weights and should be checked. Inconsistency ratio is often being used to calculate the Inconsistency in judgments, Equations (22) and (24) have been used for its calculation: where IR refers to the inconsistency ratio, I.I. is an inconsistency index, n is the number of criteria, a is the geometric mean of matrix and W (i,j) is weight vector.
Step 4: Calculation of the weighted normalized decision matrix using Equation (25): where W j represents the weight of the criteria.
Step 5: Calculation of the positive and negative ideal solution using Equations (26) and (27), respectively: where, j and J' are related to increasing and decreasing criteria, respectively, where J is associated with the positive criteria and J' is associated with the negative criteria.
Step 6: Calculation of the distance from the positive and negative ideal solution using Equations (28) and (29), respectively: Step 7: Calculation of the relative closeness to the ideal solution using Equation (30): where cl i+ is the closeness coefficient, di+ is a positive ideal solution (PIS), and d i− is negative ideal solution (NIS). The value of cl i+ ranges between 0 and 1. The larger the cl i+ value indicates the better the performance of the alternatives. In this contribution, to perform the Mathematical calculation, 500 points have been randomly selected and derived the values of GWDFs for each point, then a table was made, consisting of 16 GWDFs columns and 500 rows. Subsequently, these values have entered into SPSS and done the process. Ultimately, the TOPSIS based GWPM has been prepared considering the point values using the IDW method in the GIS environment.

Support Vector Machine (SVM)
SVM is the supervised learning system of machine learning associated with learning algorithms that analyze the data used for classification and regression analysis. It is developed by Bai et al. [59]. SVM helps in the transformations of nonlinear covariates into a higher dimensional feature space [60]. SVM is also a statistical learning theory associated with a training phase in which a training dataset of related input and target output values trains the model. The trained model will then be used to analyze a separate set of test data. SVM has two main underlying concepts for discriminating the problems, i.e., the optimum linear separating hyper-plane that separates patterns of data, and another is the kernel functions that convert the original nonlinear data pattern to a linearly separable format in a high-dimensional feature space [60]. A set of linear separable training vectors x i (i = 1, 2, . . . , n) consists of two classes, which are denoted as y i = ± 1. The SVM's goal is to find an n-dimensional hyperplane that differentiates the two groups by the total distance.
Mathematically, it can be minimized as: Subject to the following constraints: where w is the norm of the normal hyper-plane, b is a scalar base, and (w·x i ) denotes the scalar product operation. Introducing the Lagrangian multiplier, the cost function can be defined as: where λ i is the Lagrangian multiplier. It is possible to achieve the solution by double minimizing Equation (32). The standard procedures for w and b and detailed discussions can be found in Vapnik [61], Tax and Duin [62] and Yao et al. [60]. For non-separable case, one can change the constraints by setting up slack variables ξ i : Equation (32) will be modified as: where υ[0, 1] was introduced to account for misclassification [63]. Besides, a kernel function K (x i , x j ) was introduced by Vapnik [61] to account for the nonlinear decision boundary [63]. The two-class SVM method was used in this study because it was reported that Yao et al. [60] produced a more accurate map of susceptibility from the two classes of SVM. That's why Radial Basis Function (RBF) was used for kernel in this study and the two-class SVM model was first trained and then used to construct a GWPM. In this method, 1 and 0 values indicate the positive and negative relationship of groundwater occurrence. To perform the GWP mapping using the SVM, we used the ENVI 4.3. The default RBF kernel, which works well in most cases, has been used. In addition, in many studies and cases (especially in nonlinear problems), RBF provides better prediction results compared to other kernels [64]. Finally, the GWPM by SVM has been produced in the GIS.

Validation of Models
In this study, to analyze the potentiality and performance of the selected models, we have used two thresholds dependent methods i.e., ROC curve and SCAI. The ROC curve and SCAI are the significant and accurate justification methods of different models [65]. For this purpose, 30% validation and 70% training datasets have been considered by the ROC curve and SCAI methods (11). The area under curve (AUC) of the ROC method range between 0.5 to 1. If the value is nearest to 1, it indicates excellent prediction accurateness of the models [66]. The accuracy value that is AUC of the ROC is mentioned in Table 3. The AUC values have been calculated using the Equation (36). In the case of SCAI method, if the sub-class values of models decrease from very low to very high sub-classes, it indicates that models are suitable and acceptable [67]. We also used five statistical techniques in this analysis to test models' performance, including SE, SP, AC, MAE, and RMSE. Based on four possible consequences i.e., true positive (TP), false positive (FP), true negative (TN) and false negative (FN), sensitivity (Equation (37)), specificity (Equation (38)) and accuracy (Equation (39)) have been measured: TP and FP are the counts of well pixel that are correctly identified as well pixel and non-well pixel, respectively. On the other hand, TN and FN are the numbers of well pixel which are correctly classified and incorrectly classified as non-well class. SE is the ratio of the number of well pixels properly classified to the total number of well pixels predicted. SP is the ratio between the number of well pixels wrongly classified and the total non-well pixels predicted. AC is the ratio between the number of properly classified well and non-well pixels. MAE (Equation (40)) and RMSE (Equation (42)) indices have been considered to assess the disparity between the observed and predicted data. The high values of Sensitivity, Specificity, and Accuracy and low value of MAE and RMSE value indicate the good capability of the models [68][69][70][71][72]. The following five formulas have been used for statistical measures.
where X predicted and X actual is the predicted and real values in the training dataset or testing dataset of the groundwater potentiality models and n is the total number of samples in the training data set or testing dataset.

Sensitivity Analysis (SA)
It is very difficult to completely remove the uncertainty in the preparation of data layers [73][74][75]. Refsgaard et al. [75] was used different techniques e.g., Monte Carlo analysis, error propagation equations, sensitivity analysis (SA), scenario analysis, etc. for measuring the uncertainty. Sensitive analysis has been used in various studies [75][76][77] for the measurement of the effect of variable variations on model outputs, allowing then a quantitative assessment of the relative importance of uncertainty sources. In the present study, map removal sensitivity analysis (MRSA) method has been used, which was developed by Lodwick et al. [78]. The MRSA method would help to evaluate the sensitivity of the groundwater potentiality maps by removing one or more parameters from the groundwater potentiality maps. This technique has been used by several researchers to address the significant role of the effective factors [79][80][81]. It helps to identify the quantitative contribution of each groundwater conditioning factor to the uncertainty of the model output [72,82]. The percentage of contribution (PC) of each groundwater conditioning factor has been estimated by the MRSA method to explore the relative importance on the model output using the following Equation (42) [83]: where AUC all and AUC i indicate the AUC values obtained from modeling groundwater potential model using all GWDFs and the model when the ith GWDF has been excluded.

Analyzing the Multi-Collinearity (MC) of Groundwater Determining Factors
The MC problem reduces some linear models' predictive accuracy [84]. Techniques were applied in this study to assess the MC problem between GWDFs, namely tolerance (TOL) and inflation factor variance (VIF) [85]. Tolerance values of <0.1 and VIF of <10 reveal no MC problem among the GWDFs [86]. Roy and Saha [87] and Arabameri et al. [32] were used the MC test for the landslide susceptibility and groundwater potentiality mapping. The selected 16 GWDFs have been tested by SPSS. No MC problem has been found among the GWDFs, as no one value of tolerance and VIF does exceed the threshold limit (Table 4). Therefore, the selected GWDFs are suitable for the prediction of groundwater potentiality. Here, maximum tolerance and VIF values are 0.91 and 5.91 (Table 4).

Application of the Weight of Evidence (WoE)
Groundwater potentiality reclines on the positive and negative effects of the effective groundwater determining factors. The positive value of WoE indicates the chances of storage of groundwater and vice-versa. The zero value of WoE means the sub-class of factors has no role in determining the groundwater occurrence [88]. The results of WoE model have been put in Table 5. The low altitudinal zone is more potential for the accumulation of groundwater than abrupt slope and higher altitudinal areas due to the high infiltration rate and less surface runoff [89].  (Table 5). Among the five slope classes, the <2.55-degree class has the maximum value of WoE i.e., 2.97 which depicts the strong control on the occurrence of groundwater. On the other hand, the remaining four sub-classes of slope have no control over the groundwater at all ( Table 5). The north-east aspect has the highest WoE value (WoE = 1.84), which indicates the strong positive effects on the storage of groundwater. CI is the parameters of topography that reflect the elevation as a collection of convergent (channel) and divergent (ridge) areas. The CI value ranges from +100 to −100. Among the five classes, the two sub-class of CI such as <−59.21 (WoE = 2.02) and >57.64 (WoE = 2.02) has the strongest positive relationship, while others sub-layers have a negative relationship with groundwater storage (Table 5) (Table 5). Subsequently, weights have been assigned to the sub-layers of GWDFs and converted as weighted WoE layers. All weighted GWDFs have been summed up and generated a single layer of GWPM (Figure 4c). The prepared GWPM has been classified into four categories i.e., low, medium, high and very high potential zones with the help of the natural break classification method (Figure 4c). The results of GWMP by WoE model are showing that only 297.33 km 2 (15.46%) area is of very high groundwater potentiality, followed by the 583.90 km 2 (30.36%) for high, 617.37 km 2 (32.1%) for medium and 424.85 km 2 (22.09%) for low GWPZs (Table 6 and Figure 5).

Application of Random Forest (RF) Model
The RF model has been used in the present study to identify GWPZ. For carrying out the RF model, the point values based on well and non-well locations have been derived from GWDFs. The out-of-bag (OOB) of RF model is 3.32% (Table 7). The results depict that elevation (385.72), rainfall (281.63), drainage density (152.98), distance to fault (258.36), NDVI (262.47), distance to the road (189.28) and LULC (137.42) factors have great contribution in the RF process ( Figure 6). Conversely, factors such as lithology, soil type and convergence index have a tiny role in the RF process. Finally, GWPM by RF model has been prepared and classified into four classes such as low, medium, high and very GWPZs (Figure 4d). According to the GWPM, only 485.81 km 2 (25.26%) areas have very high groundwater potentiality. Other potentiality zones such as high, medium and low GWPZs are covered 9.07%, by 10.33% and 55.34% of the study area respectively (Table 6)

Application of Binary Logistic Regression (BLR)
The BLR probabilistic model has been used for GWPZ estimation. The point base data for well and non-well locations have been extracted from each GWDF. Here, BLR is expressed with binary value, i.e., 0 and 1. 1 means that the presence of well and 0 means absence of well. The coefficients of regression values have been obtained by the BLR. The results of BLR (Table 8) show that slope (0.577), soil type (6.808), lithology (2.553) and LULC (2.2942) have the reciprocal and positive impacts on the occurrence of groundwater. Among the GWCFs, elevation (0.0237), convergence index (0.0029), rainfall (0.0131), distance to fault (0.0001), distance to road (0.0002), TWI (0.2892) and TPI (0.0426) have less importance for the formation of groundwater. Conversely, Dd, distance to the river, NDVI and SPI have negative impacts on the groundwater occurrence. Afterward, the weights have been assigned to GWCFs by BLR. Finally, GWPM by BLR has been built and categorized into four categories such as low, medium, high and very high GWPZs using the natural break method. According to this classification, 15.47% of the Damghan plain has very high groundwater potentiality, followed by 38.7%, 30.91%, and 14.92% for low, moderate and high GWPZs (Figure 4a and Table 6).

Application of Binary Logistic Regression (BLR)
The BLR probabilistic model has been used for GWPZ estimation. The point base data for well and non-well locations have been extracted from each GWDF. Here, BLR is expressed with binary value, i.e., 0 and 1. 1 means that the presence of well and 0 means absence of well. The coefficients of regression values have been obtained by the BLR. The results of BLR (Table 8) show that slope (0.577), soil type (6.808), lithology (2.553) and LULC (2.2942) have the reciprocal and positive impacts on the occurrence of groundwater. Among the GWCFs, elevation (0.0237), convergence index (0.0029), rainfall (0.0131), distance to fault (0.0001), distance to road (0.0002), TWI (0.2892) and TPI (0.0426) have less importance for the formation of groundwater. Conversely, Dd, distance to the river, NDVI and SPI have negative impacts on the groundwater occurrence. Afterward, the weights have been assigned to GWCFs by BLR. Finally, GWPM by BLR has been built and categorized into four categories such as low, medium, high and very high GWPZs using the natural break method. According to this classification, 15.47% of the Damghan plain has very high groundwater potentiality, followed by 38.7%, 30.91%, and 14.92% for low, moderate and high GWPZs (Figure 4a and Table 6).

Application of TOPSIS
The TOPSIS an important MDCA approach used to delineate the GWPZs. In this research, 500 points were selected randomly and extracted point values from GWCFs. The AHP is an important knowledge-driven MDCA model, used to assign weights to GWCFs for performing the TOPSIS model. The weights of GWCFs are 0.082 (elevation), 0.088 (slope), 0.057 (aspect), 0.058 (convergence index), 0.092 (rainfall), 0.067 (drainage density), 0.063 (distance to river), 0.070 (distance to fault), 0.059 (distance to road), 0.063 (lithology), 0.050 (soil type), 0.060 (LULC), 0.061 (NDVI), 0.045 (TWI), 0.04(TPI) and 0.045 (SPI) (Figure 7). The weights of GWCFs by AHP and point base values of GWCFs have been computed using the Equations (17)- (26) and then, calculated the final weight. The GWPM by TOPSIS has been built considering the points values of GWCFs weights with the help of the inverse distance weighted (IDW) interpolation method (Figure 4b). The GWPM by TOPSIS has been classified into four classes such as low, medium, high and very high GWPZs with the help of natural break method. The results of GWMP by WoE model shows that only 290.22 km 2 (15.09%) area is very high groundwater potential, followed by the 399.46 km 2 (20.77%), 787.38 km 2 (40.94%) and 446.19 km 2 (23.2%) areas are high, medium and low groundwater potential (Table 6) Remote Sens. 2018, 10, x FOR PEER REVIEW 27 of 37

Application of TOPSIS
The TOPSIS an important MDCA approach used to delineate the GWPZs. In this research, 500 points were selected randomly and extracted point values from GWCFs. The AHP is an important knowledge-driven MDCA model, used to assign weights to GWCFs for performing the TOPSIS model. The weights of GWCFs are 0.082 (elevation), 0.088 (slope), 0.057 (aspect), 0.058 (convergence index), 0.092 (rainfall), 0.067 (drainage density), 0.063 (distance to river), 0.070 (distance to fault), 0.059 (distance to road), 0.063 (lithology), 0.050 (soil type), 0.060 (LULC), 0.061 (NDVI), 0.045 (TWI), 0.04(TPI) and 0.045 (SPI) (Figure 7). The weights of GWCFs by AHP and point base values of GWCFs have been computed using the Equations (17)- (26) and then, calculated the final weight. The GWPM by TOPSIS has been built considering the points values of GWCFs weights with the help of the inverse distance weighted (IDW) interpolation method (Figure 4b). The GWPM by TOPSIS has been classified into four classes such as low, medium, high and very high GWPZs with the help of natural break method. The results of GWMP by WoE model shows that only 290.22 km2 (15.09%) area is very high groundwater potential, followed by the 399.46 km2 (20.77%), 787.38 km2 (40.94%) and 446.19 km2 (23.2%) areas are high, medium and low groundwater potential (

Application of Support Vector Machine (SVM)
SVM is a vital machine learning data mining technique, used to recognize the potentiality of groundwater. All data layers have been reclassified into different classes using the SVM method. The SVM classification value ranges from 0 to 1, 0 indicating the absence of groundwater well in GWCFs and Conversely 1 value also indicates the presence and potentiality of groundwater formation. Among the GWDFs, the low altitude, low slopping, high rainfall, nearest distance to river, high Dd, nearest distance to road, far distance to fault, high vegetation index, H lithological units, agriculture land, arid soils, high TWI and low SPI have been considered as sub-layers of high groundwater potential, marked by the 1 values. Conversely, high altitude, slopping, low drainage density, far distance to a river, low fault distance, low vegetation index, rangeland, bare land, urban, entisols, salt flats, low TWI, and low rainfall have been identified as the 0 values because these conditions are not suitable for groundwater recharge. Thus, all GWDFs have been weighted by SVM and summed up in GIS to generate a single data layer of groundwater potential map (GWPM). The GWPM by SVM has been classified into four classes such as low, medium, high and very high GWPZs with the help of the natural break classification method (Figure 4e). The results of GWMP by WoE model show that only 360.42 km 2 (18.74%) area has very high groundwater potentiality and 248.67 km 2 (12.93%) area has low groundwater potentiality (Table 6).

Validations and Comparison of Models
Sometimes a single method of validation is not sufficient for judging the potentiality and performance of models because of the concentration of samples within a few places. Methods of AUROC, SE, SP, AC, MAE, RMSE, and SCAI were used to test the performance of models. Both training (goodness of fit) and validation (prediction accuracy) datasets have been used for judging the capability of models in producing the GWPMs of the study area. Considering the training dataset ROC curves (Figure 8  predictability to evaluate the groundwater potentiality of the Damghan sedimentary plain, although other models have good capability in mapping the groundwater potentiality.    All the statistical techniques and ROC curves used in this study for evaluating the performance of the models have judged all the models as good for mapping the groundwater potentiality in this plain. SCAI is another important validation method, used to validate the models. The SCAI values of sub-classes of all models have decreased from low potentiality to very high potentiality, indicating the appropriateness and suitability for the groundwater potentiality evaluation (Table 10). Above all, according to the threshold dependent, SCAI and statistical methods the BLR has the strongest predictability to evaluate the groundwater potentiality of the Damghan sedimentary plain, although other models have good capability in mapping the groundwater potentiality.

Sensitivity Analysis
To assess the influence of GWDFs on groundwater potentiality occurrence and to explore the effective factors with the strongest effect on the result of the groundwater potentiality prediction, a sensitivity analysis has been carried out (Table 11 and Figure 9). The results of sensitivity analysis showed in percentage contribution (PC) values of factors attained. The Pc values of the GWDFs are 7.5% (elevation), 11.35% (convergence index), 13.68% (drainage density), 11.81% (distance to road), 7.18% (distance to fault), 16.10% (distance to river), 6.19% (land use/land cover), 8.66% (lithology), 6.91% (NDVI), 9.67% (Rainfall), 7.86% (slope), 5.52% (soil), 10.10% (SPI), 12.58% (TPI), 9.51% (TWI) and 0.41% (aspect). The only slope aspect has very little contribution to the occurrence of groundwater potentiality. The results indicated that the groundwater potentiality maps of the study area are highly sensitive to elevation, lithology, drainage density, rainfall, distance to river, TPI, TWI, SPI, and distance to road. The sensitive analysis would help to reduce the variation in the model and to understand the significant geo-environmental factors that are vital for understanding the structure of model. Table 11. Sensitivity result when each factor is excluded in the binary logistic regression model.

GWDFs
Decrease of AUC (in Percentage)

Discussion
In the recent decade, the demand for water has significantly increased because of the rapid growth of population, especially in arid and semi-arid areas. The large part of Damghan sedimentary plain covering the arid and semi-arid environments groundwater is the main source of water for living. In this region, groundwater planning and sustainable management are necessary. The hydrogeologist, engineers and decision need some basic tools for managing the groundwater. GWPM may meet the basic tool of groundwater management.
GWPM is the outcome of the lithology, tectonics, topography, vegetation, rainfall, and hydrology, which are available and accessible everywhere in the environment. In this research, a different type of data has been used as the input datasets. DEMs based study provides more accurate and significant results (90)(91)(92). Different DEMs provide different results, e.g., ALOS DEM with 30 m spatial resolution provide suitable and excellent results, comparatively the ASTER and SRTM DEMs with 30 m resolution [93]. Here, the authors combined the geomorphology, geology and hydrology   Figure 9. Sensitivity result when each factor is excluded in the binary logistic regression model.

Discussion
In the recent decade, the demand for water has significantly increased because of the rapid growth of population, especially in arid and semi-arid areas. The large part of Damghan sedimentary plain covering the arid and semi-arid environments groundwater is the main source of water for living. In this region, groundwater planning and sustainable management are necessary. The hydrogeologist, engineers and decision need some basic tools for managing the groundwater. GWPM may meet the basic tool of groundwater management.
GWPM is the outcome of the lithology, tectonics, topography, vegetation, rainfall, and hydrology, which are available and accessible everywhere in the environment. In this research, a different type of data has been used as the input datasets. DEMs based study provides more accurate and significant results [90][91][92]. Different DEMs provide different results, e.g., ALOS DEM with 30 m spatial resolution provide suitable and excellent results, comparatively the ASTER and SRTM DEMs with 30 m resolution [93]. Here, the authors combined the geomorphology, geology and hydrology parameters to recognize the spatial groundwater potential. Spatial analysis is the core matter of the research for adopting the most performing approach and models for GWPMs, considering the argument topic [12][13][14][15]. Geo-environmental factors (i.e., elevation, slope, aspect, rainfall, lithology, land use/land cover, soil type, drainage density, distance to river, distance to fault, distance to road, NDVI, TWI, TPI, and SPI) were considered as the GWDFs that have been tested for the multi-collinearity problem by VIF and tolerance, and are the most effective for groundwater storage. The categorical variables such as aspect, lithology, soil type, LU/LC factors have been converted into the quantity continuous data through assigning the weight by the WoE and TOPSIS method. For the LR and RF, these GWDFs have been evaluated to prepare GWPMs taking the extracted values of GWDFs of the 500 points. The results of these models are more accurate than previous works [32]. In this work, we applied probabilistic (WoE, BLR), machine learning (SVM and RF) and multi-criteria decision approach (TOPSIS) models for building the GWPMs of Damghan sedimentary plain. These models have represented the excellent results as other works were done by Mohammady et al. [94] and Arabameri et al. [25,26]. All models, however, have very few variations in groundwater potential modeling accuracy. According to the AUROC, SE, SP, Accuracy, MAE and MRSE among the five models, the BLR models (for training dataset AUC = 0.933, SE = 0.852, SP = 0.828, AC = 0.839, MAE = 0.216 and RMSE = 0.314 and for validation dataset AUC = 0.943, SE = 0.909, SP = 0.846, AC = 0.875, MAE = 0.214 and RMSE = 0.311) have better capability for mapping the groundwater potentiality than other models. Recognizing the significance of each variable for groundwater storage is very difficult. The soil, lithology, altitude, rainfall, LU/LC, NDVI, Dd, distance to fault factors are dominant factors among 16 GWDFs for the formation of groundwater. The SA is depicting the contribution in producing the uncertainty in the GWPM and the factor distance from the river has the highest contribution to the variation of output of model. Similar to the Shahroud plain, the Damghan sedimentary plain regions consists of the large bare land, rangeland, and urban land, interrupting the water infiltration into sub-surface layer, while agriculture land with aquifer locations are receiving the larger water into the sub-surface and also signify the hydrologic properties [95]. According to TOPSIS model, the rainfall, slope, elevation, LU/LC, soil type factors have been highly prioritized by the AHP model, suggesting the most potential for groundwater formation. Such findings are confirmed with the work of Arabameri et al. [32]. Among the 16 GWDFs, the elevation is the most important topographic component that influences the groundwater recharge. In fact, at a lower segment, the Damghan sedimentary plain is almost flat, where water stagnation and associated infiltration of water is maximized. On the contrary, high altitudes, associated with open and v-shaped slopes promote runoff due to local physiography. The methods applied for validating the GWPMs are showing outstanding accuracy, and ensemble models have better capabilities than the individual modes. Such ensemble models have been shown to be more reliable in this analysis than the other models used by the researchers for GWPM in various other locations [96,97]. The proper methods can have the ability to produce GWPMs, and that can be used for planning purposes. The used probabilistic, machine learning and ensemble models have excellent accuracy and may be used for groundwater management in this plain region.

Conclusions
Today, GWPM is an effective groundwater resource management method. Through the Over-extraction of groundwater in the low groundwater, the potential region can be limited by the GWPM. With the advancement in the technical field, different techniques for the spatial modeling of groundwater are introducing day by day. So, it very difficult to say what method would be best for spatial modeling. However, in the present research, five methods (BLR, TOPSIS, WoE, RF, and SVM) have been used for modeling the groundwater and the compared among them to answer the question of what model is relatively better for the Damghan sedimentary plain. The GWPM approaches are more appropriate to predict the potential of groundwater. The GWPMs have been produced with the help of RS and GIS techniques. RS and GIS both combinedly helped to perform the works such as identification of well, thematic data generation, classification, and final map generations. The RS and GIS-based study are cost and time saving, accurate, and provide meaningful results. The R studio is an important machine learning program that helps to perform different kinds of models such as logistic regression, random forest, naive Bayes tree, support vector machine, artificial neural network, and several other methods. R program based model performance is more easy, accurate, efficiency and perfect. The GWPMs, produced by the selected methods have been categorized into four categories i.e., low, medium, high and very high potential classes. The results of GWPMs show that the very high potentiality zones are covered with by an area of 290.22 km 2 (15.09% by TOPSIS), 297.34 km 2 (15.46% by WoE), 485 km 2 (25.26% by RF), 297.53 km 2 (15.47% by BLR) and 360.42 km 2 (18.74% by SVM) out of 1923.27 km 2 areas. The worthiness of GWPMs has been significantly validated by the ROC and SCAI methods and five statistical measures i.e., SE, SP, AC, MAE, and MRSE. According to the results of the ROC, SCAI methods and statistical measures, these models are excellent for the prediction of GWPZ. The very high or excellent GWPZs have been found in the low elevated and less sloppy area. The arid soils are covered by high potentiality of groundwater. In the case of the LU/LC and vegetation index, the agriculture land and high vegetation density areas have high potentiality of groundwater. Conversely, high altitude, sloppy land, urban area, rangeland, salt flats, entisols soil type have the low potentiality of groundwater formation. The GWPMs may be used as tool in this study area for managing and developing the groundwater. The resulting maps can also assist decision-makers, planners, and engineers in choosing the ideal location, groundwater distribution for further groundwater exploration. Therefore, Damghan sedimentary plain region has high potentiality of groundwater storage which can be saved by sustainable use, obstructing groundwater pollution, increasing the people's awareness and suitable government policy regarding the amount and way water use.