Groundwater-Potential Mapping Using a Self-Learning Bayesian Network Model: A Comparison among Metaheuristic Algorithms

Owing to the reduction of surface-water resources and frequent droughts, the exploitation of groundwater resources has faced critical challenges. For optimal management of these valuable resources, careful studies of groundwater potential status are essential. The main goal of this study was to determine the optimal network structure of a Bayesian network (BayesNet) machine-learning model using three metaheuristic optimization algorithms—a genetic algorithm (GA), a simulated annealing (SA) algorithm, and a Tabu search (TS) algorithm—to prepare groundwater-potential maps. The methodology was applied to the town of Baghmalek in the Khuzestan province of Iran. For modeling, the location of 187 springs in the study area and 13 parameters (altitude, slope angle, slope aspect, plan curvature, profile curvature, topography wetness index (TWI), distance to river, distance to fault, drainage density, rainfall, land use/cover, lithology, and soil) affecting the potential of groundwater were provided. In addition, the statistical method of certainty factor (CF) was utilized to determine the input weight of the hybrid models. The results of the OneR technique showed that the parameters of altitude, lithology, and drainage density were more important for the potential of groundwater compared to the other parameters. The results of groundwater-potential mapping (GPM) employing the receiver operating characteristic (ROC) area under the curve (AUC) showed an estimation accuracy of 0.830, 0.818, 0.810, and 0.792, for the BayesNet-GA, BayesNet-SA, BayesNet-TS, and BayesNet models, respectively. The BayesNet-GA model improved the GPM estimation accuracy of the BayesNet-SA (4.6% and 7.5%) and BayesNet-TS (21.8% and 17.5%) models with respect to the root mean square error (RMSE) and mean absolute error (MAE), respectively. Based on metric indices, the GA provides a higher capability than the SA and TS algorithms for optimizing the BayesNet model in determining the GPM.

rithms [34,35] or group techniques [36][37][38] have recently been used. The Bayesian network (BayesNet) model has been successfully used in hazard modeling so far [39,40], but one of the disadvantages of this model is the inability to determine the optimal network structure.
The main objective of this research is to develop GPM using a BayesNet model combined with genetic algorithm (GA), simulated annealing (SA), and Tabu search (TS) algorithms to determine the optimal network structure of the BayesNet in the Baghmalek Plain, Iran. Second, the study aims to compare the prediction capabilities of the different metaheuristic algorithms. For this purpose, the GA, SA, and TS algorithms have been used to optimize the network structure of the BayesNet model to prepare a groundwater-potential mapping. To our knowledge, the BayesNet-GA, BayesNet-SA, and BayesNet-TS models have not been investigated to address this issue yet.

Materials and Methods
This research was accomplished in five levels. In the first level, a spatial database was created, consisting of the locations of the springs and the criteria affecting the potential of groundwater. In the second level, the important criteria affecting the potential of groundwater was determined using the OneR technique. In the third level, the spatial relationship between the effective criteria and locations of the springs was achieved using the certainty factor (CF) method. In the fourth level, spatial modeling was performed using ensembles of the BayesNet model and metaheuristic (i.e., GA, SA, and TS) algorithms. In the last level, the GPM was evaluated using metric indicators. A flowchart of our research is shown in Figure 1.

Study Area
The town of Baghmalek is in the east of Khuzestan Province, Iran. It is geographically located at 49 • 53 E and 31 • 31 N, at an altitude of 917 m above sea level. The temperature in the center of the study area in summer reaches a high of 42 • C, and in winter reaches a low of 0 • C. The average annual rainfall is 500 mm. Due to the high use of groundwater in this area for agriculture, in recent years there has been a clear decrease in groundwater levels. The location of the study area, along with the distribution of the springs, is shown in Figure 2. The town of Baghmalek is in the east of Khuzestan Province, Iran. It is geographically located at 49° 53' E and 31° 31´ N, at an altitude of 917 m above sea level. The temperature in the center of the study area in summer reaches a high of 42°C, and in winter reaches a low of 0 °C. The average annual rainfall is 500 mm. Due to the high use of groundwater in this area for agriculture, in recent years there has been a clear decrease in groundwater levels. The location of the study area, along with the distribution of the springs, is shown in Figure 2.

Spring Inventory Map
Our general assumption was that the presence of groundwater in regions with high groundwater yields was highly likely [26]. A total of 187 springs were chosen, with a yield greater than 11 m 3 /h and a mean pH and electrical conductivity (EC) of 7.1 and 490 μmhos/cm, respectively. In this study, the locations of the springs were provided by the Water Resources Management Organization (WRMO) of Iran. Google Earth was used to verify the locations of the springs. Of the locations, 70% (131 locations) were randomly used for modeling (i.e., training), and 30% (56 locations) were randomly used for the modeling evaluation (i.e., validation).

Criteria Affecting the Potential of Groundwater
The topography, meteorology, hydrology, human activities, and environmental features from the earth are major criteria affecting the potential of groundwater among the different elements [13,22]. To determine the GPM, 13 factors were selected in our research. Considering the previous articles, these criteria included: altitude, slope angle, slope aspect, plan curvature, profile curvature, topography wetness index (TWI), distance to river, distance to fault, drainage density, rainfall, land use/cover, lithology, and soil [26,35]. ArcGIS 10.3 and SAGA GIS software were used to prepare and classify the effective criteria. The following is a description of each of the effective criteria. The method for preparing and classifying the effective criteria is summarized in Appendix Table A1.
Altitude: Because of differences in soil and vegetation caused by altitude differences, groundwater is usually found at low altitudes. The shuttle radar topography mission (STRM) digital elevation model (DEM) was downloaded from https://earthex-

Spring Inventory Map
Our general assumption was that the presence of groundwater in regions with high groundwater yields was highly likely [26]. A total of 187 springs were chosen, with a yield greater than 11 m 3 /h and a mean pH and electrical conductivity (EC) of 7.1 and 490 µmhos/cm, respectively. In this study, the locations of the springs were provided by the Water Resources Management Organization (WRMO) of Iran. Google Earth was used to verify the locations of the springs. Of the locations, 70% (131 locations) were randomly used for modeling (i.e., training), and 30% (56 locations) were randomly used for the modeling evaluation (i.e., validation).

Criteria Affecting the Potential of Groundwater
The topography, meteorology, hydrology, human activities, and environmental features from the earth are major criteria affecting the potential of groundwater among the different elements [13,22]. To determine the GPM, 13 factors were selected in our research. Considering the previous articles, these criteria included: altitude, slope angle, slope aspect, plan curvature, profile curvature, topography wetness index (TWI), distance to river, distance to fault, drainage density, rainfall, land use/cover, lithology, and soil [26,35]. ArcGIS 10.3 and SAGA GIS software were used to prepare and classify the effective criteria. The following is a description of each of the effective criteria. The method for preparing and classifying the effective criteria is summarized in Appendix A Table A1.
Altitude: Because of differences in soil and vegetation caused by altitude differences, groundwater is usually found at low altitudes. The shuttle radar topography mission (STRM) digital elevation model (DEM) was downloaded from https://earthexplorer.usgs. gov (accessed on 28 February 2021) at a spatial resolution of 30 m × 30 m. Topographic parameters were obtained from the DEM in ArcGIS10.3 and SAGA GIS 2.1.2 software, including altitude, profile curvature, slope aspect, slope angle, and plan curvature (Figure 3a) [41].
Slope angle: The slope angle affects vertical water infiltration and surface runoff velocity. In areas with lower slopes, runoff recharge is also higher (Figure 3b) [42].
Slope aspect: The amount of water in the soil is affected by the slope aspect, and the physiographic trends and the direction of precipitation are related to the slope aspect ( Figure 3c) [22].

CF Method
To solve for uncertainty and heterogeneity of the input data for modeling, the CF method was implemented. This method was developed by Shortliffe and Buchanan [49]. The weight of each category of criteria that affected the potential of groundwater was calculated using Equation (2) [49]. TWI: The TWI shows the location and size of saturated source areas, and the production of surface runoff is affected by topographic/morphological conditions. The TWI is calculated using Equation (1) (Figure 3d) [43].
where A s . is the specific surface of the catchment and α is the slope gradient. Plan and profile curvature: This controls the speed of water flow. The values of the profile curvature represent the morphology of the topography, with a positive curvature attesting to an upward concave curvature, a negative curvature attesting to an upward convex curvature, and a value of zero attesting to a flat surface. For a longer period of rainfall, a concave slope retains more water (Figure 3e,f) [26].
Distance to river: This is the primary source of groundwater recharge in semiarid regions that affects the potential of groundwater (Figure 3g) [42].
Distance to fault: The movements of groundwater springs over an area can be controlled by various faults (Figure 3h) [44].
Drainage density: The drainage characteristics of each basin control the subsurface hydrological status of each area. The drainage density varies in different regions due to land use effects, different climatic regimes, and natural landscapes (Figure 3i) [45].
Rainfall: One of the most essential and natural sources of groundwater recharge that impacts the infiltration for the potential of groundwater is the amount of rainfall. A 30-year average rainfall dataset for 22 meteorological stations obtained by the WRMO was used to prepare the rainfall map. The rainfall map was constructed in ArcGIS 10.3 by using the kriging interpolation method (Figure 3j) [17].
Lithology: Permeability and porosity are two criteria that affect the movement and occurrence of groundwater as its efficiencies increase with growing permeability and porosity. Lithology plays a crucial role in permeability and porosity ( Figure 3k) [46].
Soil: Surface-water permeability is controlled by soil properties and is directly related to permeability (Figure 3l) [47].
Land use/cover: Groundwater movements are affected by diversity in land use/cover. Evapotranspiration, surface runoff, and groundwater recharge are affected by land use/cover. For the preparation of the land use/cover map, Landsat-7 images from 2013 and 2019 were used. To this purpose, 420 training points obtained via a Global Positioning System (GPS) were used to create a land use/cover map, of which 70% were used for training, and the remaining 30% were used to test the accuracy of the map. The land use/cover maps were produced at an accuracy of 90% using the maximum likelihood algorithm in ENVI 4.8 software (Figure 3m) [48].

CF Method
To solve for uncertainty and heterogeneity of the input data for modeling, the CF method was implemented. This method was developed by Shortliffe and Buchanan [49]. The weight of each category of criteria that affected the potential of groundwater was calculated using Equation (2) [49].
where cf is the weight of each category, pp a is the conditional probability of the existing number of springs in a category, and pp s is the initial probability of the existence of all springs in study area. The favorability values (pp a , pp s ) were obtained from overlaying the existing spring distribution layer in GIS with each data layer and measuring the frequency of occurrence of spring. The value of pp a is given by the ratio of the spring area falling into a certain category and the total area of that category. The value of pp s can be determined by dividing the total area of the spring with the total area of the study. The value in the CF method is between −1 and +1. Positive values indicate an increase in certainty, and negative values indicate a decrease in certainty. In this research, the CF method was used in order to match and take into account the effect of different classes of independent parameters in the modeling [49].

Baysian Network (BayesNet) Machine-Learning Model
The BayesNet machine-learning model is a knowledge representation tool that enables effective management of dependency/independence relationships between the random variables that make up the modeling-problem domain [50]. This model consists of two components, structure and parameter. The structure of the BayesNet machine-learning model, of nodes that represent system variables and E ⊂ → V × → V is a collection of edges. It describes the relationship of direct dependence between these variables [39]. For each the set of parents of the variable x i . The cumulative distribution on → V can be calculated using Equation (3) [18,50].

Genetic Algorithm (GA)
A GA that mimics biological evolution can be used to solve optimization problems with and without restrictions. This algorithm was introduced by Holland in 1975 as one of the first metaheuristic algorithms [51,52]. In this algorithm, each individual solution is called a chromosome, and consists of single arrays called genes. This algorithm repeatedly modifies a population of solutions. At each step, the GA randomly selects individual solutions from the present population and uses them as parent chromosomes to create offspring chromosomes for the new generation. Genes are transformed from current chromosomes into new generation chromosomes using operators. After creating successive generations, the population moves toward an optimal solution [53,54].

Simulated Annealing (SA) Algorithm
The SA algorithm is a simple and effective metaheuristic search used to solve hybrid optimization problems [55]. The SA algorithm simulates a gradual refrigeration process to solve an optimization problem. The objective function of such a problem is similar to the energy of a substance, which must be minimized by defining a virtual temperature. The SA algorithm is a probabilistic algorithm in which a solution is proposed to exit the local optimization. Also, it is a memoryless algorithm in the sense that there is no mechanism for storing information during the search [56,57].

Tabu Search (TS) Algorithm
The TS algorithm is a direct search algorithm for optimizing complex nonlinear problems. In this algorithm, the transition from a possible initial solution to a possible secondary solution is called a move [58]. The list of the last points to be examined is called the taboo list, which is proposed based on human memory. The purpose of addressed list is to prevent relocation to areas that have already been surveyed. This can expand the algorithm search area. First, the algorithm randomly selects an initial answer, then all neighboring points are identified, and the value of objective function is calculated [59].

Hybrid Model
One of the main components in a BayesNet machine-learning model is the structure of the network. The traditional methods for determining the best network structure in their search space can acquire local optimization and ignore global optimization [60]. For this purpose, metaheuristic algorithms are an effective solution to determine the best network structure. In our research, to optimize the structure of the BayesNet machinelearning model, metaheuristic (i.e., GA, SA, and TS) algorithms were used. The objective function of optimizing the BayesNet structure (B s ) for the D dataset using metaheuristic algorithms is defined in Equation (4): where P(B s ) is the previous network structure, Γ(.) is the gamma function, and N ij and N ijk are Dirichlet distribution indices for a set of parameters. The hybrid model was applied using the Waikato Environment for Knowledge Analysis (WEKA 3) software. To implement the models, the data was prepared in ArcGIS 10.3 software and then transferred to the WEKA 3 software for modeling. Finally, the output of WEKA 3 software was transferred to ArcGIS 10.3, and the final groundwater potential map was prepared (Figure 4).

OneR Technique
In this study, the OneR technique was used to determine the feature selection and select the most effective criterion in modeling the groundwater-potential map [61]. The OneR technique is a one-tier decision tree that contains a set of rules. It is simple and often provides good rules for characterizing data structures. It also uses the computational error

OneR Technique
In this study, the OneR technique was used to determine the feature selection and select the most effective criterion in modeling the groundwater-potential map [61]. The OneR technique is a one-tier decision tree that contains a set of rules. It is simple and often provides good rules for characterizing data structures. It also uses the computational error ratio and several rules to gain the weight of each effective criterion [62].

Validation Indices
To evaluate the hybrid models, Precision, Recall, Specificity, Sensitivity, Kappa, receiver operating characteristic (ROC)-area under the curve (AUC), root mean square error (RMSE), and mean absolute error (MAE) indices were used. The ROC curve is a graph in which the occurrence or nonoccurrence of springs correctly predicted by the model is plotted on the horizontal axis (1-Specificity) versus incorrectly predicted pixels (Sensitivity) on the vertical axis. The area below the ROC curve is called AUC. Its value varies between 0.5 and 1; the closer it is to one, the higher the modeling efficiency is. These indicators were calculated using Equations (5)-(11) [61,63].
where TP is the number of spring values that were correctly predicted, FP is the number of spring values that were incorrectly detected, TN is the number of spring values that were rejected correctly, FN is the number of spring values that were incorrectly rejected, Y is the real value,Y is the predicted value, and N is the number of samples. The RMSE and MAE indices were used to measure the error between the actual and predicted values [64,65].

Determining the Important Criteria
The importance we determined for the criteria that affect the potential of groundwater is shown in Figure 5. Based on the results of OneR technique, the criteria of altitude, lithology, and drainage density had the most significant impact on the occurrence of groundwater in study area, which is consistent with previous articles [66,67]. Water 2021, 13, x FOR PEER REVIEW 12 of 23 Figure 5. The importance of the effective criteria determined using the OneR technique.

Results of the CF Method
The spatial relationship using the CF method between the criteria affecting the potential of groundwater and the location of springs is presented in Table 1. Based on the results of altitude criterion, the class of 724-1084 m provided the highest correlation with the presence of groundwater (0.355). As the altitude increased, the potential of groundwater decreased. Also, the reason for this phenomena can be expressed as the reduction of aquifer thickness at high altitudes [26]. Considering the results of slope angle criterion, the slope angle class of less than 8° supplied the highest weight (0.220). The potential of groundwater decreased with an increasing slope. The main reason for this behavior is that the nature of water tends to descend from points at high altitudes to lower altitudes, and also to move in the direction of the slope and accumulate in the lowlands [14]. In the results for the slope aspect criterion, the south direction with a weight of 0.392 supported the highest importance. Due to receiving more rainfall and humidity than other areas, it played a more significant role in controlling the formation of groundwater [41]. Regarding the results of plan curvature criterion, the concave class yielded the highest weight and the most significant impact on the occurrence of groundwater (CF = 0.427) [23]. Recognizing the results of profile curvature criterion, the class of less than −0.0054 produced the highest weight, equal to 0.608. With an increasing TWI, the speed of water passage from high altitudes will decrease, and water will have more opportunity to penetrate the ground. Therefore, the potential of groundwater will increase [68]. But in this study, the middle classes of TWI had a greater impact on the occurrence of groundwater. One of the reasons for the negative weight of the category larger than 6.1 is that a small number of springs (5) are located in this category. The results show that values between 4 and 6 for this index are more important in the occurrence of groundwater, which is consistent with [4]. For the results of distance to river, the class of 0-100 m yielded the highest weight (CF = 0.744) [69]. For distance to fault, the class of 0-300 m reported the highest value (CF = 0.474). Faults can transfer groundwater in hard formations. The presence of a fault facilitates infiltration conditions. Therefore, at distances closer to the fault, the potential of groundwater is higher [14]. In the results for the drainage density criterion, the maximum weight was related to the class of 0.487-0.731, with a value of 0.618. As the drainage density increased, the area's permeability increased, and thus provided suitable conditions for groundwater recharge [16]. For the results of rainfall criterion, the class of 550-500 mm had the highest weight (CF = 0.490). With increasing rainfall, the probability of groundwater occurring increased. Since the study area has a semiarid climate, the average rainfall

Results of the CF Method
The spatial relationship using the CF method between the criteria affecting the potential of groundwater and the location of springs is presented in Table 1. Based on the results of altitude criterion, the class of 724-1084 m provided the highest correlation with the presence of groundwater (0.355). As the altitude increased, the potential of groundwater decreased. Also, the reason for this phenomena can be expressed as the reduction of aquifer thickness at high altitudes [26]. Considering the results of slope angle criterion, the slope angle class of less than 8 • supplied the highest weight (0.220). The potential of groundwater decreased with an increasing slope. The main reason for this behavior is that the nature of water tends to descend from points at high altitudes to lower altitudes, and also to move in the direction of the slope and accumulate in the lowlands [14]. In the results for the slope aspect criterion, the south direction with a weight of 0.392 supported the highest importance. Due to receiving more rainfall and humidity than other areas, it played a more significant role in controlling the formation of groundwater [41]. Regarding the results of plan curvature criterion, the concave class yielded the highest weight and the most significant impact on the occurrence of groundwater (CF = 0.427) [23]. Recognizing the results of profile curvature criterion, the class of less than −0.0054 produced the highest weight, equal to 0.608. With an increasing TWI, the speed of water passage from high altitudes will decrease, and water will have more opportunity to penetrate the ground. Therefore, the potential of groundwater will increase [68]. But in this study, the middle classes of TWI had a greater impact on the occurrence of groundwater. One of the reasons for the negative weight of the category larger than 6.1 is that a small number of springs (5) are located in this category. The results show that values between 4 and 6 for this index are more important in the occurrence of groundwater, which is consistent with [4]. For the results of distance to river, the class of 0-100 m yielded the highest weight (CF = 0.744) [69]. For distance to fault, the class of 0-300 m reported the highest value (CF = 0.474). Faults can transfer groundwater in hard formations. The presence of a fault facilitates infiltration conditions. Therefore, at distances closer to the fault, the potential of groundwater is higher [14]. In the results for the drainage density criterion, the maximum weight was related to the class of 0.487-0.731, with a value of 0.618. As the drainage density increased, the area's permeability increased, and thus provided suitable conditions for groundwater recharge [16]. For the results of rainfall criterion, the class of 550-500 mm had the highest weight (CF = 0.490). With increasing rainfall, the probability of groundwater occurring increased. Since the study area has a semiarid climate, the average rainfall of 500 mm had a more significant impact on the occurrence of groundwater. In the study area, more rainfall occurred in areas with a high altitude and slope, which are not prone to groundwater occurrence; therefore, the rainfall criterion was affected by these two criteria, and high amounts of rainfall (>550 mm) had negative weight. Regarding the results for the lithology criterion, the MuPlaj unit showed the highest weight (CF = 0.639). The MuPlaj unit, owing to its high permeability, had the most significant impact on the potential of groundwater in the study area [14]. Contemplating the results of the land use/cover criterion, the agricultural class provided the highest weight (CF = 0.347) and the most significant impact on the occurrence of groundwater. Agricultural irrigation water, mainly supplied from river water, penetrates the saturated aquifer and feeds the aquifer during irrigation [70]. Finally, for the soil criterion, the Inceptisols class suggested the highest weight (CF = 0.370). The reason for this is that the increased permeability of Inceptisols provides the necessary conditions for the formation of aquifers in areas covered with this type of soil [26].

Results of the Hybrid Models
To optimize the BayesNet machine-learning model, a spatial database provided the points of presence of spring (value 1) and the absence of spring (value 0). The weights were obtained from all criteria by the CF method. Nonoccurrence points were randomly created as equal to the training and validation points. The control parameters in the metaheuristic algorithms are presented in Table 2. The diagrams for the comparison of training and validation data with target data (occurrence and nonoccurrence of springs) are shown in Figure A1. The evaluation results of the hybrid models using metric indicators for the training and validation data are summarized in Table 3. As observed in Table 3

GPM Using Hybrid Models
After training (calibrating) the hybrid models, to generalize the hybrid models for the entire study area, the entire area first was converted into points, and the values produced by the hybrid models were prepared for the entire study area. We used ArcGIS 10.3 software to prepare the GPM using hybrid models. The GPM was divided into five categories: very low potential, low potential, medium potential, high potential, and very high potential, using the natural breaks classification method. The GPM using four models is shown in Figure 6. According to the BayesNet-GA model of GPM, the percentages for the very low, low, medium, high, and very high potential categories were 2%, 24%, 37%, 27%, and 10%, respectively. The GPM using the BayesNet-SA model showed that the percentages for the very low, low, medium, high, and very high potential categories respectively were 18%, 21%, 24%, 24%, and 13%, respectively. In the BayesNet-TS model, the percentages for the very low, low, medium, high, and very high potential categories were 15%, 24%, 25%, 24%, and 12%, respectively. The GPM using the BayesNet model showed percentages for the very low, low, medium, high, and very high potential categories of 20%, 21%, 25%, 20%, and 14%, respectively.

Validation of GPM
To evaluate the results using the ROC-AUC, 30% of the occurrence data and the data of the absence were used for the individual springs. The results of the GPM evaluation with three hybrid models using the ROC-AUC are shown in Figure 7 and Table 4. According to the results, the ROC-AUC's values for the BayesNet-GA, BayesNet-SA, BayesNet-TS, and BayesNet models were 0.830, 0.818, 0.810, and 0.792, respectively. The results showed that the BayesNet-GA model was more accurate than the other three models.

Validation of GPM
To evaluate the results using the ROC-AUC, 30% of the occurrence data and the data of the absence were used for the individual springs. The results of the GPM evaluation with Water 2021, 13, 658 15 of 20 three hybrid models using the ROC-AUC are shown in Figure 7 and Table 4. According to the results, the ROC-AUC's values for the BayesNet-GA, BayesNet-SA, BayesNet-TS, and BayesNet models were 0.830, 0.818, 0.810, and 0.792, respectively. The results showed that the BayesNet-GA model was more accurate than the other three models.   The results of the three hybrid models indicated a superior accuracy in preparing the GPM. Therefore, it can be concluded that metaheuristic algorithms (i.e., GA, SA, and TS) were very useful in optimizing BayesNet model to prepare GPM. Among the metaheuristic algorithms, the GA was more accurate than the other two algorithms. The use of objective-function values to perform the optimization process without the need for additional information such as function derivatives, the simplicity of the search process, and the flexibility can be demonstrated as the best advantages of the GA application [71]. After the GA, the SA algorithm had the next-highest accuracy in optimizing the BayesNet model. The advantage of the SA algorithm was very low memory consumption and the ability to pass local optimization due to a randomized process [72]. However, this algorithm was less accurate in optimizing the BayesNet model than the GA due to its high dependence on the initial values of parameters and the lack of proper determination for the initial temperature parameter [73]. The TS algorithm may stop at a local optimization. Still, it is not clear whether the result obtained is a local optimization or a global optimization, and the local optimization obtained may depend on the initial response [74]. Therefore, the TS algorithm was less accurate in optimizing the BayesNet model than the other two algorithms.  The results of the three hybrid models indicated a superior accuracy in preparing the GPM. Therefore, it can be concluded that metaheuristic algorithms (i.e., GA, SA, and TS) were very useful in optimizing BayesNet model to prepare GPM. Among the metaheuristic algorithms, the GA was more accurate than the other two algorithms. The use of objective-function values to perform the optimization process without the need for additional information such as function derivatives, the simplicity of the search process, and the flexibility can be demonstrated as the best advantages of the GA application [71]. After the GA, the SA algorithm had the next-highest accuracy in optimizing the BayesNet model. The advantage of the SA algorithm was very low memory consumption and the ability to pass local optimization due to a randomized process [72]. However, this algorithm was less accurate in optimizing the BayesNet model than the GA due to its high dependence on the initial values of parameters and the lack of proper determination for the initial temperature parameter [73]. The TS algorithm may stop at a local optimization. Still, it is not clear whether the result obtained is a local optimization or a global optimization, and the local optimization obtained may depend on the initial response [74]. Therefore, the TS algorithm was less accurate in optimizing the BayesNet model than the other two algorithms.

Conclusions
Groundwater plays a vital role in human activities in various sectors, such as agriculture and industry. Therefore, determining the potential of groundwater using a machinelearning model is an effective step in groundwater management. To achieve higher accuracy in determining the potential of groundwater and solving the BayesNet model's problem, metaheuristic algorithms were used. Based on the results, metaheuristic algorithms improved groundwater-potential mapping with the BayesNet model. Based on the results of metric indices and the ROC-AUC values for the GPM, the GA had a higher accuracy than the SA and TS algorithms. According to the RMSE and MAE criteria, the BayesNet-GA model improved the GPM estimation accuracy by 4.6% and 7.5% (21.8% and 17.5%) compared to the BayesNet-SA and BayesNet-TS models. The results of the CF method showed that the criteria of distance to river and lithology had a higher weight in modeling compared to the other criteria. Also, the criteria of altitude, lithology, and drainage density were of higher importance in groundwater in study area. Prepared maps can help manage and control groundwater resources in arid and semiarid regions.

Conflicts of Interest:
The authors declare no conflict of interest.