Novel Ensemble of MCDM-Artificial Intelligence Techniques for Groundwater-Potential Mapping in Arid and Semi-Arid Regions (Iran)

: The aim of this research is to introduce a novel ensemble approach using Vise Kriterijumska Optimizacija I Kompromisno Resenje (VIKOR), frequency ratio (FR), and random forest (RF) models for groundwater-potential mapping (GWPM) in Bastam watershed, Iran. This region suffers from freshwater shortages and the identification of new groundwater sites is a critical need. Remote sensing and geographic information system (GIS) were used to reduce time and financial costs of rapid assessment of groundwater resources. Seventeen physiographical, hydrological, and geological groundwater conditioning factors (GWCFs) were derived from a spatial geo-database. Groundwater data were gathered in field surveys and well-yield data were acquired from the Iranian Department of Water Resources Management for 89 locations with high yield potential values ≥ 11 m 3 h − 1 . These data were mapped in a GIS. From these locations, 62 (70%) were randomly selected to be used for model training, and the remaining 27 (30%) were used for validation of the model. The relative weights of the GWCFs were determined with an RF model. For GWPM, 220 randomly selected points in the study area and their final weights were determined with the VIKOR model. A groundwater potential map was created by interpolating the values at these points using Kriging in GIS. Finally, the area under receiver operating characteristic (AUROC) curve was plotted for the groundwater potential map. The success rate curve (SRC) was computed for the training dataset, and the prediction rate curve (PRC) was calculated for the validation dataset. Results of RF analysis show that land use and land cover, lithology, and elevation are the most significant determinants of groundwater occurrence. The validation results show that the ensemble model had excellent prediction performance (PRC = 0.934) and goodness-of-fit (SRC = 0.925) and reasonably high classification accuracy. The results of this study could aid management of groundwater resources and assist planners and decision makers in groundwater-investment planning to achieve sustainability.

Many approaches have been used to identify areas with high groundwater potential and the array of approaches has increased in recent years because of the growing pressure to clarify the distribution of groundwater and to identify new sources of water. High-precision maps are needed to better manage aquifers that are in use and to avoid costly exploration for new deposits, especially during droughts. Hybrid groundwater mapping techniques are a vital necessity [35]. Ensemble machine learning frameworks can be used to develop better prediction models and they easily cope with complex multi-dimensional data [40]. Kordestani et al. [30] hybridized an EBF and BRT (i.e., creating an EBF-BRT) algorithm and demonstrated a groundwater contamination-susceptibility model. Naghibi et al. [47] hybridized four data mining models and FR into a novel ensemble model for a groundwater-related study. Their ensemble model exhibited better performance than standalone applications of each separate model. Chen et al. [35] developed a novel hybrid approach with WoE, LR, and functional tree (FT) models that individually outperformed the models. Miraki et al. [40] used a random subspace (RS) ensemble-RF (i.e., RS-RF) classifier to model groundwater susceptibility and found that the hybrid model outperformed individual models. Razavi-Termeh et al. [48] integrated three bivariate statistical models with RF and logistic tree models to map groundwater potential and revealed that the RF statistical models improved prediction accuracy. Though the idiographic studies demonstrated hybrid models that performed strong groundwater susceptibility modeling in given study areas, there is no evidence to show that they are universally appropriate or applicable. The applicability of the novel approaches for better prediction of groundwater susceptibility models needs further evaluation.
Therefore, the main purpose of this research is to evaluate the integration of MCDM with machine-learning approaches to map groundwater potential. FR models, based on simple statistical methods, are used to examine correlations between groundwater-related factors and groundwater potential. This ensemble approach is novel and has not yet been applied in groundwater potential mapping. This study evaluated the performance of the Vise Kriterijumska Optimizacija I Kompromisno Resenje (VIKOR)-RF-FR ensemble approach for groundwater potential mapping. Statistical evaluation metrics-area under the receiver operating characteristic curve (AUC-ROC), FR, and seed cell area index (SCAI)-were used to determine the goodness-of-fit and predictive performance of the ensemble algorithm. Data from the advanced land-observing satellite (ALOS) phased array type L-band synthetic aperture radar (PALSAR) digital elevation model (DEM) with a spatial resolution of 12.5 m were used to acquire topographical and hydrological factors. ALOS-PALSAR DEM has been determined to have a higher accuracy than advanced spaceborne thermal emission and reflection radiometer (ASTER) and shuttle radar topography mission (SRTM) DEMs [49].
MCDM models are useful for solving complex problems because of their capacity to enable examination and selection of alternatives based on criteria [50][51][52][53][54][55][56][57][58][59]. MCDMs can analyze basic statistical trends in big data in various study areas [60,61]. These models are an important part of modern decision-making science which is decision-making in the face of varying decision criteria and selection from alternative outcomes. MCDM methods have been developed not only because of the diversity of processes and problems that require assessment of multiple factors, but because experts desire to improve decision making with better techniques [62][63][64][65]. The improvements are the outcomes of advances in mathematical science, computational science, and computer technology. Machine-learning models have also gained popularity among environmental managers and decision makers [66][67][68]. Machine-learning models compare favorably to other models because they do not require a defined and strict set of a priori assumptions as do many other statistical methods. These models are also used to assess the relative importance of each of the conditioning factors [66]. The FR model simplifies the probabilistic relationships between the dependent and independent variables to quantify the relationships between the reclassified factors and groundwater capacity. The relationships calculated using the FR model enabled easier interpretation of groundwater potential indicated by the rating of each factor.

Study Area
Bastam watershed, with an area of 1329.43 km 2 , is located between 36°27′02″ and 36°47′13″N and between 54°24′23″ and 55°11′08″E in Semnan Province, Iran ( Figure 1). The region is topographically mountainous with an elevation range from 1357 to 3893 m.a.s.l. and an average elevation of 2078.3 m asl. The average slope in the study area is 13.5° and the maximum is 70.6°. The region's climates are arid and semi-arid: the average annual rainfall of the study area is 262.4 mm and the average annual temperature is 12.8 °C [69]. Valley terrace deposits, stream channels, braided channels, flood plain deposits, low elevation piedmont fans, and valley terrace deposits are the most important lithological units of the study area [70]. Entisols (32.9%), entisols/inceptisols (32.4%), and mollisols (21.5%) are the most widespread soils in the study area [71].
Geologically, Devonian-Eocene rocks, including andesitic and basaltic volcanic, conglomerate, and sandstone, are distributed as mountain ranges in the northern part of the study area extending from the eastern to the western boundaries. Cretaceous granitic and volcanic rocks are intruded as mountain ranges and isolated hills. They are found in the eastern part of the study area. Moreover, Quaternary sediments, fluvial deposits along streams, and collovium, are found in the central section of study area. Groundwater in this area is stored in these Quaternary sediments, both alluvial and sedimentary formations, in cracks, fractures, faults, bedding planes, and large cavities of consolidated rocks of Permo-Carboniferous meta-sediment and granitic aquifers.
Due to Iran's aridity, industrialization, and population growth, groundwater is a vital natural resource. More than 90% of groundwater is used for agriculture [72]. During the period 1971-2014, the number of water wells in the country increased dramatically, from 47,000 to more than 789,000. There are currently about 196,000 deep wells and 593,000 semi-deep wells in Iran and they extract an average of 0.17 and 0.02 million m 3 , respectively [72]. Groundwater consumption in Iran has gone from less than 18 billion m 3 in 1971 to more than 61 billion m 3 in 2014. In recent years, average groundwater levels in the country have dropped by about 51 cm, a decrease of 5.086 billion m 3 /year [48]. Approximately 20,000 people depend on the groundwater resources of the Bastam watershed for personal consumption and agriculture [73]. Mismanagement of water in this region has led to a 15 m drop of water tables over a 20-year period and to land-surface subsidence [73].

Methodology
This research involved six steps ( Figure 2): (i) Prepare the database-A groundwater-inventory map (GWIM) was created and 17 groundwater conditioning factors (GWCFs) were added as layers to a GIS (ArcGIS). The groundwater inventory locations were separated into two randomly selected groups for training and validation.
(ii) Assess the data-A multi-collinearity analysis was performed on the GWCFs using indices of tolerance (TOL) and the variance inflation factor (VIF). TOL < 0.2 and VIF > 5 indicate a multicollinearity issue among the variables [74]. The linearity of the GWCFs was assessed, and increasing linearity between the factors will reduce the accuracy of the models [65].
(iii) Calculate relative importance of GWCFs-Because a VIKOR MCDM model cannot determine the relative importance of predictive variables, a RF model was used. To do this, the same number of non-well locations were randomly selected in the study area as there were wells yielding more than 11 m 3 /h and they were used to train the RF model. Next, the Extract Multi Values to Point tool in ArcGIS extracted the values of each GWCF at well (coded 1) and non-well (coded 0) locations. The data were formatted in Excel and imported into R 3.1.1 software. The Rattle package in R software was used to calculate the weight of each GWCF. Given a groundwater dataset, GWD = (xi, yi), (i = 1, 2, 3, ... N) with xi ∈ R M are the conditioning factors. N is the total number of samples, and R M is the total number of conditioning factors, yi ∈ (1,0) is the output which contains two classes, well and nonwell locations, with the first coded "1" and the second coded "0". The bootstrap aggregating algorithm generated n subsets, which consist of m factors (m ≤ M). A CART algorithm constructed a tree classifier for each subset. The tree classifiers were aggregated to an RF classifier. Two parameters of each the RF, n, and m, were determined when building the RF model. The parameter, n, should be large enough to ensure the robustness of the RF model (i.e., n = 500).
(iv) Prepare the decision matrix-Two hundred and twenty-two random points ( Figure 3) (called alternatives) were selected in the study area and their GWCF values were extracted. A decision matrix of 17 columns (the GWCFs) and 220 rows (alternatives) was created. FR models were used to calculate the values of the alternatives' pixels. In MCDM procedures, GWCF values should have either ascending or descending orders. The actual values did not follow this trend in this study, so FR was used to fit an ascending trend. Higher GWCF values indicated locations that were more likely to have groundwater.
(v) Run the hybrid model and map groundwater-The VIKOR was run to develop the groundwater potential model. The produced values indicated strength of association with the presence of groundwater at the 220 locations. Using the Ordinary Kriging tool in ArcGIS 10.5, the distribution of groundwater potential was mapped. Ordinary Kriging has been shown to efficiently map nonlinear relationships and has been used for a variety of real-world applications [75].
(vi) Assess performance-The receiver operating characteristics (ROC) curves and their area under the curve (AUROC), frequency ratio (FR), and seed cell area index (SCAI), were used to evaluate the accuracy of the model predictions.

Data Preparation
The most important step in groundwater-potential mapping (GWPM) is the creation of a wellinventory map (WIM) [64]. Extensive field surveys were conducted, and well-yield data were acquired from the Iranian Department of Water Resources Management to characterize the known spatial distribution of groundwater in the region. Actual pumping-test production measurements of wells expressed in m 3 h −1 were used to represent yield. Eighty-nine wells with a high yield (≥ 11 m 3 h −1 ) were selected for GWPM [1] and modeling ( Figure 1). The data from these wells were randomly grouped into a training dataset of 62 (70%) wells and a validation dataset of 27 (30%) to be used to test the model [76].
Geo-environmental features were used as GWCFs to predict groundwater potential [1]. Based on a literature review [1,30,32,36,37,39,40], availability of data, environmental characteristics, scale of research, and a multicollinearity test, 17 RS-derived GWCFs that reflected topography, hydrology, and LULC were selected. Some GWCF measurements were derived from the ALOS PALSAR DEM with a spatial resolution of 12.5 m [1,52,53]. RS and thematic maps were also used to measure others. The GWCFs are elevation, slope, aspect, stream power index (SPI), topographic position index (TPI), topographic wetness index (TWI), terrain ruggedness index (TRI), convergence index (CI), distance to streams, distance to faults, distance to roads, drainage density (DD), rainfall, soil type, land use and land cover (LULC), lithology, and normalized difference vegetation index (NDVI) (Figure 4a-q).
Elevation, slope, and aspect were derived from ALOS PALSAR DEM in ArcGIS. Distances to streams, roads, and faults were computed as Euclidean distances pixel centroids to the nearest polyline. The vector-based information was then rasterized to coincide with the resolution of the DEM. Road and fault-factor maps (scales of 1:50,000 and 1:100,000, respectively) were obtained from the National Geographic Organization of Iran (www.ngo-org.ir) and the Geological Society of Iran (GSI) (http://www.gsi.ir/), respectively. Descriptions of lithological units in the study area were compiled. A factor map of streams was extracted from the DEM in ArcHydrov10.4 and ArcGISv10.4. Soil texture, lithology, and LULC were derived from three thematic maps: the soil map of the Semnan Agricultural and Natural Resources, Research and Education Center, the geological map of the Geological Society of Iran (GSI) (http://www.gsi.ir/), and the land use map created by the Iranian Soil Conservation and Watershed Management Research Institute (https://www.scwmri.ac.ir). All three of these maps were at the 1:100,000 scale.
Annual precipitation data from nine weather stations (at Mojen, Farahzad, Bastam, Abr, Karkhaneh, Semnan, Shahroud, Tarzeh, and Shah Kuh-e Bala) for the period from 1986 to 2016 were mapped and Kriging [1,24,51,54] was used to map annual rainfall patterns in ArcGIS 10.5 [69]. Land uses significantly influence the distribution of groundwater [1]. Greater vegetative cover tends to both indicate and promote the presence and extent of aquifers, as forests and denser vegetation reduce runoff rates and encourage infiltration more than do barren and sparsely vegetated lands. Equations (1)-(4) were used to calculate TWI, SPI, TRI, and TPI [77][78][79].
where A S is the specific catchment area of the basin (m 2 /m), β is slope steepness in degrees, x indicate elevation of each neighbor cell to cell (0, 0) (m), and max and min show the largest and smallest elevation value among nine neighbor pixels, respectively. Epixel depicts the elevation of the cell, and Esurrounding is the mean elevation of the neighboring pixels. TPI compares the elevation of each pixel in the DEM with the average height of the pixels around it. This factor enables classification of landscapes into morphological classes. The positive and negative values indicate that a pixel is higher or lower than the pixels that surround it [80]. The plan curvature reflects the directional variations along a curve. The effect of plan curvature on slope erosion is that it causes convergence and divergence of water along and away from a flowing stream [81]. The profile curvature represents the amount of elevation variation along the flow path. When the concavity of the surface curvature is increasing, its values are negative, and, if decreasing, they are positive [82]. This index indicates the flow velocity, erosion (in negative amounts), sedimentation (in positive values), and geomorphology of the area [83]. TPI, plan curvature, profile curvature, and convergence index were extracted from the DEM in system for automated geoscientific analyses (SAGA-GIS).
NDVI is used to evaluate vegetation types and density using quantitative analysis of albedo or reflectance of visible and NIR portions of the electromagnetic spectrum. On a scale from −1 to +1, high vegetation density usually has values between 0.05 and 1, bare soil areas are close to 0, and water bodies have values closer to −1. Vegetation depends upon moisture and promotes infiltration into soils and therefore is highly indicative of a likely presence of groundwater [1]. NDVI values were computed from a LANDSAT-8 image with 30 m spatial resolution acquired from the USGS (https://earthexplorer.usgs.gov/) ( Table 1).
The GWCFs used in the GWPM were converted to raster layers at a 12.5 m resolution. The raster matrix was 2395 columns by 1460 rows.
Classes of lithology: A) limestone, B) marl and shale, C) tuffaceous shale, Andesitic and basaltic volcanics, D) limestone, and shale, E) marl, shale and detritic limestone, F) Red conglomerate and sandstone, G) siltstone and shale, P) fluvial conglomerate, and sandstone, I) braided channel and flood plain deposits, J) dolomite, shale and sandstone.

Frequency Ratio (FR)
FR is the probability of occurrence of a specific phenomenon. FR determines the relationship between GWCFs and well locations [61]. To calculate FR for a conditioning factor, Equation (5) is used [84]: where A is the number of wells in a class, B is the total number of wells, C is the number of pixels in that class, and D is the number of pixels [85].

VIKOR (Vlse Kriterijumsk Optimizacija Kompromisno Resenje)
The VIKOR method was first introduced by Opricovic [86]. VIKOR is derived from compromising programming. This model is an MCDM method that determines the best solution for a discrete decision that has disproportionate and conflicting criteria. This method focuses on ranking and selection from alternatives and it determines a compromise solution to guide decision-making. The compromise solution is a feasible solution that is the closest solution to the ideal solution. Compatibility is also based on concessions [87]. The following six steps are involved in the VIKOR model: (1) Prepare the decision matrix.
(2) Calculate the normalized matrix, as shown in Equation (6): where, xij is the ith alternative performance of the jth criterion, and m is the alternative numbers.
(3) Calculate a weighted normalized matrix, as shown in Equation (7): where Wj is the weight of the criterion. (4) Identify the ideal positive (Equation (8)) and negative (Equation (9)) options: where f i + is the positive ideal solution for the ith criterion and f i − is the negative ideal solution for the ith criterion. (5) Calculate the utility index (Equation (10)) and incompatibility index (the distance from positive and negative ideal solution) (Equation (11)): where S i and R i are the distance of the ith alternative to the ideal positive and negative solutions, respectively.
(6) Calculate the VIKOR index and determine the final weight of the alternatives, as shown in Equation (12): where V is a constant (0.5), S * is min Si, S − is max Si, R * is min Ri, and R * is max Ri. The greater the value of the VIKOR index in one alternative, the greater the importance of that alternative.

Random Forest
RF is one of the most well-known machine-learning algorithms [88]. RF is nonparametric and based on decision trees. Many decision trees are grown in the classification of RF algorithms [89]. The Expands each tree completely without pruning.
The GWCFs are prioritized using the mean decrease accuracy and the mean decrease Gini. Using the mean decrease accuracy, as opposed to the mean decrease Gini, determines the priority of effective factors more effectively, especially when environmental factors are being compared [68].

Validation of Results
AUC-ROC, FR, and SCAI criteria help to discern the quality of definite and probable identification, as well as the quality of prediction [65]. AUC-ROC represents a model's predictive accuracy by describing its ability to predict whether predefined events have occurred or will occur [90]. AUC-ROC curves indicate the sensitivity of the model to the percentage of cells or unstable units correctly predicted by the model against the percentage of predicted unstable cells within the total. These values reflect the model's ability to correctly distinguish between positive and negative observations in a validation sample. High sensitivity indicates that most predictions were correct (true positives), while specificity values of 1 indicate a high number of false positives. In the AUC-ROC, the false positive rate (1 specificity) is shown on the x axis (Equation (13)) and the true positive rate (sensitivity) on the y axis (Equation (14)): where TN is true-negative, FP is false-positive, TP is true-positive, and FN is false-negative [85]. The quantitative-qualitative relationship classes between AUC-ROC and prediction accuracy, ranging from 0 to 1, are: excellent (0.9-1), very good (0.8-0.9), good (0.7-0.8), moderate (0.6-0.7), and weak (0.5-0.6) [80]. FR and SCAI criteria were used to evaluate the separations between classes. These criteria determine the classification accuracy of a model. The FR is the ratio of the percentage of wells in each class to the percentage area of that class [91], and SCAI is the ratio of the percentage area of each zoning class to the percentage of wells in each class [92].

Multi-Collinearity Analysis
The multicollinearity test shows that there is no collinearity among conditioning factors ( Table  2). TOL among factors ranges from 0.15 to 0.84 and VIF ranges from 1.1 to 4.8. Therefore, all GWCFs were included in the GWPM process. * VIF is variance inflation factor.

Determining the Relative Importance of GWCFs using the RF Model
The RF model with an out-of-bag error rate = 8.57% was applied in R using the Caret package. This scenario implies that the accuracy of the model is equal to 91.43%, and exhibits the results of the confusion matrix. According to the results (Table 3), of the 80 non-well locations, 74 (92.5%) were predicted to be non-well and 6 (7.5%) were predicted to be locations of wells. Of the 80 well locations, 11 (13.75%) were predicted to be non-well locations and 69 (86.25%) were predicted to be wells. The determination of the relative importance of the GWCFs (using mean decrease accuracy of the RF model) ( Figure 5) shows that LULC (18.5), lithology (13.9), and elevation (11.9) were the main factors that identified the locations of greater amounts of groundwater. These were followed by TPI (9.1), NDVI (7.9), TWI (7.7), distance to road (7.4), drainage density (6.3), slope aspect (5.7), soil type (5.3), distance to fault (5.03), TRI (4.9), slope degree (4.8), rainfall (4.7), CI (4.6), distance to stream (2.4), and SPI (1).

Figure 5.
Relative importance of conditioning factors using random forest model.

Determining the Value of Each Pixel of GWCFs using the FR Model
The FR model revealed the spatial relationships between the GWCFs and wells ( Table 4). Analysis of the values of each GWCF at each well location (Figure 6a-q) shows that locations with lower elevations, slope angles, CI, SPI, distance to road, distance to stream, rainfall, TRI, and NDVI had lower FR values (i.e., these factors have strong inverse relationships with groundwater potential). Locations with higher TWI and drainage density have higher values of FR (i.e., strong positive relationships to groundwater potential). In terms of LULC and lithology, agriculture and quaternary formations, not surprisingly, were strongly associated with well locations.

Application of VIKOR Model
Extracting the GWCF values at random locations, the constructed decision matrix, and the final weight of random locations using the VIKOR model indicate that groundwater potential is strongly indicated by weighting values. The weights of random locations range from 0 to 1. Higher values indicate greater groundwater potential.

Validation of Results
Validation of the ensemble model is based on the prediction rate curve (PRC) and the success rate curve (SRC) (Figure 9). The ensemble model, with PRC = 0.93 and SRC = 0.92, has excellent groundwater potential prediction accuracy. The FR and SCAI model results ( Figure 10) show ascending and descending trends, respectively. In terms of statistical assessment of the modeling values, FR increased and SCAI decreased, indicating that the classification accuracy of the ensemble model was reasonable.

Discussion
Groundwater resources are crucial in regions with limited surface water, particularly in arid and semi-arid climates [81]. When these regions have significant amounts of industrial and agricultural water demands, aquifers can be consumed at rates that exceed recharge rates and are at increasing risk of irreparable damage. Therefore, it is important to have accurate assessments of the spatial dimensions of groundwater in a watershed. Such assessments can improve water-management planning and use strategies [93]. Groundwater mapping can be improved and made more economical by identifying the geophysical and hydrological conditioning factors that are associated with subsurface storage. Several methods have been developed to map groundwater potential, and each has advantages and disadvantages.
This research developed a novel state-of-the-art approach using an ensemble of statistical, machine-learning, and MCDM approaches in conjunction with RS and GIS techniques to map groundwater potential. MCDM models, statistical models, and machine-learning algorithms (like VIKOR, FR, and RF, respectively) can be used to support scientific decisions for problems that have a diversity of management criteria [94]. VIKOR can analyze both quantitative and qualitative data, uses simple mathematical formulas, produces stable results, is flexible, is highly efficient, is very accurate, it can factor uncertainty into analyses, is comprehensive in its analysis, produces visual representations of the data, enables comparison of alternatives, and can incorporate and analyze large amounts of data. FR was used for its simplicity of equations, easily interpretable results, and its high efficiency [95]. RS and GIS save time and are cost-effective ways to improve data gathering and analysis [96][97][98][99].
Machine-learning algorithms can model natural phenomena (e.g., groundwater potential and groundwater presence) with nonlinear relationships. They do not require removal of outliers, data transformations, or statistical assumptions. They can handle complex nonlinear relationships between conditioning factors and groundwater potentiality and can automatically analyze effectivefactor (i.e., predictor) interactions. RF was the machine-learning model employed in this research because of its advantages. The RF method can accommodate different types of predictor variables as well as missing values. It can assess the relationships and interactions between effective factors. It can handle uncertainty, can apply several input factors without eliminating any, and can return a small set of categories to maintain high prediction accuracy [54,67]. Classification accuracy, however, is affected by the number, scale, types, and precision of inputs. Incorporating all suitable factors increases modeling accuracy. Moreover, compared to other models, RF can include more datasets [53]. It can overfit some datasets with noisy classification and regression tasks. It can be biased toward attribute data with more categories, and it can also be biased toward attributes with more levels. Therefore, the variable-importance scores from RF are not reliable for such data. There are a few drawbacks to this algorithm, however. It does require significant computational memory, analysis of large datasets consumes a lot of time, there are temporal high costs for pruning, there are a high number of overlapping end nodes, and errors accumulate at each layer as trees grow.
The analysis of the spatial relationships between GWCFs and well locations using the FR model showed that elevation, slope, CI, SPI, distances to the nearest road, distances to the nearest stream, less rainfall, TRI, and NDVI are inversely correlated to groundwater potential. TWI and drainage density are directly correlated to well locations. These relationships are consistent with previous studies [29,33]. Elevation and relief play important roles in groundwater potential because weather and climatic conditions vary greatly between elevations and cause differential soil development and vegetation growth patterns. Areas at low elevations, with gentle slopes, and with low drainage density are usually associated with highly permeable material, higher vegetative cover, and less relief. These areas have higher groundwater potential. Conversely, higher elevations, steeper slopes, and higher drainage density usually have less soil, more impermeable subsurface material, sparse vegetation, and runoff-promoting relief. TWI, CI, TRI, and SPI topographic indices predict soil moisture, groundwater flow, and slope stability. These indices have been used to reflect spatial patterns of soil-moisture [77]. In terms of LULC and lithology, agricultural areas and regions with Quaternary formations have the strongest relationships to wells. Agriculture tends not to occur in places without access to water and wells are not often drilled in places that lack demand for water. The LULC categories covary with soil conditions and subsequently reflect groundwater occurrence [100]. Lithology influences groundwater because lithological variation is associated with varying porosity and hydraulic conductivity in rock formations and subsequent soils. Soil has an impact on recharge due to permeability factors and composition, and so it can influence the potential size of an aquifer [100].
The RF model indicated that LULC, lithology, and elevation were most strongly related to groundwater flow rates. This is consistent with other research [1,8,31,34,38,100]. Rahmati et al. [8] mapped groundwater potential in the Mehran Region of Iran by combining RF and maximum entropy models using data from 163 wells. They analyzed ten GWCFs and elevation, drainage density, lithology, and land use were found to be the greatest determinants of groundwater production. Arabameri et al. [1] combined RF, WoE, binary logistic regression, and the technique for order-preference by similarity to ideal solution (TOPSIS) for groundwater potential mapping on the Shahroud Plain in Iran. Analysis of the relationships between 122 wells and 15 GWCFs revealed that LULC, soil type, and slope were key correlates to groundwater. Naghibi et al. [38] used boosted regression tree, classification and regression tree, and RF machine-learning models for GWPM in Iran. Thirteen hydrological-geological-physiographical conditioning factors were assessed at 864 springs. Elevation and drainage density were the most important correlates to spring locations.
Comparing the goodness-of-fit and performance of the ensemble model using AUROC, FR, and SCAI methods revealed that the ensemble model had excellent prediction accuracy, and this was consistent with other research [30,43]. Combined or ensemble models are often more effective than are individual, stand-alone models. Kordestani et al. [30] showed that by combining EBF and BRT for GWPM, prediction accuracy (AUROC) climbed from 75% for EBF alone to 82% for EBF-BRT. Four data-mining models and FR were combined as an ensemble model [43] and showed better performance by the ensemble model because it reduced overfitting.
Overfitting occurs due to noise, limited training data, and classification complexity. Solutions are difficult because of the width of the array of potential causes. Three aspects of this study suggest that overfitting could have occurred: the area from which the non-well locations were selected was confined to the region that contained wells even though the study region extended well beyond the region covered by wells, the ratio of 70:30 training to validation datasets was used without an assessment of the accuracy of the sampling ratio, and despite the multicollinearity assessment, GWCFs still had noise.
A lack of confidence (uncertainty) might result from doubt about the completeness of knowledge. Data could be suspected of being unclear, inaccurate, unreliable, inconclusive, or potentially false [101][102][103][104]. Machine-learning techniques like RF can account for nonlinear relationships and can handle uncertainty in data. In an RF model, the prediction interval has two main components that determine its width: one represents the uncertainty of the predicted mean (or some other parameter)-this is the confidence interval-and the other represents the variability of individual observations around that mean. BRT uses uncertainty reduction to arrive at a final number of trees [105]. There are many sources of error in datasets used for modeling: measurement errors, sampling bias, limitations in field data collection, genetic variability, etc. These errors can affect model accuracy [106]. The strength of the novel ensemble model introduced here to model natural phenomena with nonlinear relationships is confirmed by the results of this study.
There are several things that can enhance GWPM accuracy and they should be the foci of future research. First, there should be a complete analysis of the implications of limiting or extending the area from which the selection of non-well locations is made to places where wells are located or to areas beyond the distribution of wells. Second, the conventional 70:30 ratio of the datasets used for training and validation should be revisited, as proportion may be causing issues. Third, the methods used to select features may be having unseen effects on the modeling process, and fourth, there should be a comparison of the VIKOR-RF-FR ensemble to other ensemble models to predict groundwater locations.

Conclusions
Many methods have been used to assess regional groundwater potential. This study created an MCDM statistical machine-learning ensemble model to spatially predict groundwater in Bastam Watershed, Iran. VIKOR, FR, and RF models were combined to map and assess the characteristics of groundwater-bearing landscapes. Like many other parts of the world, the Bastam Watershed suffers from a shortage of water. Seventeen conditioning factors (elevation, slope, aspect, SPI, TPI, TWI, TRI, CI, distance to stream, distance to fault, distance to road, drainage density, rainfall, soil type, LU/LC, lithology, and NDVI) were assessed at 89 well locations to determine each conditioning factors' statistical association with groundwater presence and productivity. RF modeling revealed that LULC, lithology, and elevation were the most important factors predicting groundwater potential and production. Validation confirmed that the ensemble model had high prediction accuracy. The VIKOR-RF-FR ensemble model is a systematic, objective, and scientific decision-making model that can be used to support groundwater resource development and management plans. A GWPM can help planners guide future groundwater exploration. The same analysis in other districts with similar topographic and geological conditions can save time and money in their quests for protecting and developing groundwater resources. A limitation of this research is that a fixed combination of GWCFs were considered for modeling. Different combinations of GWCFs should be explored in order to further improve groundwater potential modeling.