Reciprocal analysis of groundwater potentiality and vulnerability modeling in the Bahabad Plain, Iran

Heuristic and statistical groundwater quality assessment models are efficient tools in the zoning of groundwater vulnerability to contamination. An important reciprocal methodology, yet neglected in Iran, was conducted to assess the performance of three groundwater vulnerability models, namely GODS, SI, and DRASTIC, and a data mining model for groundwater potential, maximum entropy (MaxEnt). For both the training and validation stages for the MaxEnt model, the Mahalanobis distance technique was adopted. The vulnerability rates obtained from the DRASTIC model with a coefficient of determination (R2) value of 0.76 had a statistically significant correlation with nitrate concentrations in the 21 wells, compared to SI and GODS. The DRASTIC model can better reflect the vulnerability of groundwater resources to contamination. The impact of the vadose zone with an average effective weight of 33 is more important than other parameters, followed by depth than groundwater (D) (32.01), net recharge (R) (28.95), and the aquifer media (A) (18.1). These weights may not be changed. MaxEnt showed significant performance in both the training and validation stages with the respective area under the receiver operating characteristic curve (AUROC) values of 0.907 and 0.901. A reciprocal analysis between the vulnerability map in the superior model and the groundwater potential map derived from MaxEnt revealed that areas with high groundwater potential are still in safe state but require more attention as the top priority for amendment practices. In addition, approximately 8.7% of the entire study area has a high vulnerability to contamination, which requires immediate pragmatic actions.


Introduction
Nowadays, many countries face the problem of water scarcity (and therefore, groundwater is one of the most important sources that provides essential water for drinking and ecosystem (Zhao and Pei 2012;Singh et al. 2015;Qian et al. 2020). At least two billion people depend on groundwater as the main source of their drinking water (Li and Merchant 2013).
Knowledge and understanding of the impacts of groundwater withdrawal on long-term trends in water quality are essential for the effective and sustainable management of groundwater resources (Qian et al. 2020). Groundwater vulnerability modeling is critical for implementing programs to protect groundwater quality (Li and Merchant 2013). However, various man-made activities are often responsible for deteriorating the quality of groundwater systems (Datta et al. 2011;Gnanachandrasamy et al. 2020). Furthermore, characteristics of the type of bedrock in a particular way affect the quality of water resources (Salamati et al. 2017 modeling of groundwater vulnerability enables managers to achieve better risk-based decision-making and consequently increases the reliability and credibility of groundwater vulnerability predictions (Rahmati et al. 2019). On the one hand, due to the susceptibility of groundwater resources and the difficulties involved in designating protection measures against contamination, they have drawn due attention from the scientific community (Arzu Firat and Fatma 2013). The groundwater quality in a region is regarded as a function of natural and anthropogenic factors . Therefore, water resource managers are concerned about the effects of contaminants, e.g., impacts of nitrate levels on human health and ecological processes (Li and Merchant 2013).
Generally, the quality of groundwater resources depends on the recharge sources of the aquifer, precipitation, and geochemical processes. It is also affected by massive contaminations charged from anthropogenic activities (Reza and Singh 2010;Srinivas et al. 2013). However, in arid regions due to low precipitation, the risk of groundwater resource contaminations imposed by anthropogenic agents is more discernible. Therefore, it is imperative to monitor the quality of valuable groundwater resources to find practical solutions to estimate and improve the water quality (Hashmi et al. 2009;Adimalla and Venkatayogi 2017).
Among pollutants, nitrate has severe impacts on human health and the environment throughout hydrological systems (Panagopoulos et al. 2011;Rahmati et al. 2019;Stefania et al. 2019). Aquifer contamination by contaminates is also severely affected by topographical attributes (Jahanshahi et al. 2020). On the bright side, many models have been developed to tackle such data and give good results on the groundwater contamination status.
Due to the complexity of monitoring groundwater contamination, the traditional methods face difficulties in extracting useful information about contamination sources Zhang et al. 2020). The methodologies developed incorporating numerical simulation models also encounter implicit constraints (Datta et al. 2011). However, it is inevitable to use numerical models to assess and prevent imminent groundwater contamination (Takizawa 2008;Locatelli et al. 2019;Rahmati et al. 2019). During recent decades, a variety of models/technics have been developed for modeling groundwater contamination. There are several overlapping index methods, such as WQI (Ghosh and Goel 2017), DRASTIC (Boumaiza et al. 2021), GODS (Ghazavi and Ebrahimi 2015), and SI (Ghouili et al. 2021) for zoning groundwater resource quality. Several studies have been conducted on the importance of examining the quality of groundwater resources in arid/semi-arid areas using different index models (Hosseini and Saremi 2018;Ng'anga et al. 2019;Bahir et al. 2020). Therefore, this study aimed to investigate the applicability and efficiency of the GODS, SI, and DRASTIC models on the interpretation of groundwater quality and zonation of aquifer contamination in one of the most important groundwater resources in the Bahabad Plain, in central Iran.
Additionally, a one-sided study, based solely on groundwater contamination status without knowing the actual spatial potential of groundwater, leads to partial and, sometimes, misleading messages. Most studies, especially those conducted in Iran, are based on methodologies that can merely address aquifer vulnerability to contamination; however, reciprocal groundwater potentiality and contamination modeling (GPCM) is often neglected. Such an analysis can answer whether areas with high vulnerability to contamination also have a high groundwater potential or not, and therefore different patterns can entirely change management plans and decision-making strategies. Depending on data availability, many models have been proposed to pinpoint the areas with high groundwater potential. Data mining models can extract the emerging pattern stored in different geoenvironmental data and prove helpful in groundwater potential mapping (Rahmati et al. 2016;Razandi et al. 2017;and Raney and Leopold 2018). Due to the nature of groundwater resources and drawing on the literature review, the maximum entropy was used to map the groundwater potential map of the Bahabad Plain and complete the missing piece of our interactive analysis. MaxEnt was coupled with the Mahalanobis distance (MD) method to have an enriched dataset for the training and validation stages. The vulnerability of groundwater to contamination as regional threats and the spatial potential of groundwater as an opportunity would give us a holistic view on where to be concerned and where to protect.
Based on the premises mentioned above, the main objectives of this study include the following: (1) modelling the vulnerability of aquifers to contamination using DRASTIC, SI, and GODS models, (2) modelling groundwater potentiality using the coupled MaxEnt-MD model, and (3) integrating the latter two modelling stages into the GPCM framework to attain a holistic blueprint for decision making and management.

Study area
The study was carried out in the Bahabad basin, Central Iran (Fig. 1). The area is between 55° 42′ and 56° 17′ N latitudes and 31° 32′ and 32° 13′ E longitudes. The study area extends for about 2162 km 2 with respective minimum and maximum elevations of 1083 and 3052 m.a.s.l. The site is characterized by an arid climate with an average annual precipitation of 81.5 mm and an average annual temperature of 17.9 °C. In the last decade, the precipitation amount is decreased to 68 mm, which indicates an alarming condition.
The poor-condition rangeland with about 53% areal coverage is the primary land cover of the study area, followed by rock outcrops (34.6%), bare lands (6.8%), and agriculture (3.1%) while the remaining area is shared between moderate-condition rangelands, orchards, residential, and salty lands, accounting for 2.3% of the study area (www. yzrw. ir). The main water usage in the area corresponds to agriculture, water-intensive industries, and drinking, which are mainly delivered by either public suppliers (domestic deliveries) or private self-supplied sources (mostly wells). The annual water supply from the installed wells amounts to 42.19 MCUM (million cubic meters), out of which 82% alone is demanded by agriculture, followed by 12% for industry, and 6% for drinking. Therefore, a reciprocal assessment of groundwater potentiality and vulnerability is essential. The methodological flowchart of this study is presented in Fig. 2.
In terms of geology, Bahabad Plain is a part of central Iran and the Lut desert basin. The existing geological strata in Bahabad Plain are very diverse, with different facies. Figure 3 depicts the geological map of the plain. In terms of stratigraphy, it is diverse, and various facies have been formed, from the oldest Precambrian sediments to the most recent alluvial sediments and desert areas of the Quaternary period. Volcanic and plutonic masses with different facies are exposed in the western heights of the region. From the viewpoint of expansion of different geological facies, the most significant expansion of rock outcrops (27.40%) belongs to Paleozoic deposits. Paleogene and Neogene sediments with a 25.98% expansion, as well as Jurassic shale sandstone with a 24.35% expansion, have created half of the stone outcrops in the region (Fig. 3).

Data collection
All the required parameters of physical attributes of the primary aquifer as the main groundwater resource in the region were obtained from all the 21 operating wells in the area. The parameters are explained in detail in the following sections, obtained from Iran Water Resource Management Company. All required data examined in 21 wells are interpolated using the inverse distance weighting (IDW) method in ArcGIS 10.8.

DRASTIC
Although a number of spatial models designed to assess groundwater pollution hazard have been proposed, DRASTIC is arguably the model most widely used for such efforts. Aller et al. (1985) offer a detailed account of the DRASTIC methodology, its evolution, and guidelines for applications. The model was designed to be a simple, easyto-use, nationally applicable tool for groundwater pollution hazard.
For aquifer vulnerability assessment, the DRASTIC model incorporates seven factors, including depth to groundwater (D), net feed or net recharge (R), middle aquifer environment (aquifer media) (A), soil media (S), topography (slope) (T), the impact of the vadose zone (unsaturated zone) (I), and hydraulic conductivity (C), following Eq. 1.
where DI stands for DRASTIC Index while w and r are the rank and weight of each factor, respectively (Kazakis and Voudouris 2015). In this method, the vertical vulnerability of groundwater resources is based on the assumption that (1) DI = D r D w + R r R w + A r A w + S r S w + T r T w + I r I w + C r C w Fig. 2 Methodological flowchart adopted in this study. RSP slope position, SPI stream power index, TPI topgrophic position, TRI terrain ruggedness index, TWI topographic wetness index contamination occurs on the ground surface and then percolates into the saturated zone concurrently with rainfall. In addition, contamination travels with water at the same rate.

GODS
The GODS method for determining groundwater vulnerability is basically divided into three general categories namely process method based on simulation, statistical method, and overlapping method.
The GODS method is a simple, practical, and experimental approach that gives a quick assessment of the contamination sources (Foster 1987;Huang et al. 2013). The GODS calculates the vulnerability index of aquifers in response to contaminations using four indices, including the type of aquifer (G), lithology of the unsaturated zone (O), water table depth (D), and soil properties (S), following Eq. 2 (Paez 1990;Shrestha et al. 2017).

SI
Where SI denotes the Sensitivity Index and r and w indices are the rank and weight of the layer, respectively. Nitrate concentration in 21 wells was acquired from Iran Water Management Company and validated vulnerability maps derived from the three models mentioned above. In where x o is the measured NO 3 value, xp is the predicted value (obtained from the DRASTIC, GODS, and SI models), and n is the number of samples. x 0 represent the average of the measured NO 3 . The value of R 2 is always between zero and one. The closer it is to one, the stronger the correlation between the two datasets.

Weight assignment and relative importance of the input parameters of the vulnerability model
The results of weight assignment and relative importance of the input parameters of the DRASTIC and SI models are presented in Table 1. According to Foster (1987) and Paez (1990), equal weights were assigned to GODS model parameters. The assigned weights to DRASTIC model parameters are based on the values suggested by Aller et al. (1987). Furthermore, weights were assigned to SI parameters based on the values suggested by Stigter et al. (2006) and Ribeiro et al. (2017). For instance, according to Aller et al. (1987) on applying the DRASTIC model to 10 different counties across the United States, the parameter of the impact of Vadose zone (I) with an average effective weight of 33 is relatively more important than other parameters, followed by the depth to groundwater parameter (D) with an average effective weight of 32.01, the net recharge parameter (R) with an average effective weight of 28.95, and aquifer media (A) with an average effective weight of 18.1. It is recommended that these weights may not be changed. (4) Groundwater potential mapping using the coupled MaxEnt-MD Separate modelling of groundwater potential is required to grasp a holistic view of the degree to which contamination sources can threaten groundwater sources in the study area.
In other words, since the primary aquifer in the region is of an unconfined type, sites other than those subjected to well logging are expected to have groundwater potential. Moreover, if we rely merely on wells to conceive the groundwater potential, it would best give us a two-class presence-absence groundwater potential map that distorts the reality of nature; that is, the region exhibits different degrees of potentiality. Therefore, a probabilistic approach would be required. In this regard, data mining tools can extract the emerging pattern of interest, handling different types of data, and the validity of their results can be tested by the pieces of evidence observed in nature. Data mining models are mainly categorized into presenceonly and presence-absence; both offer magnificent outputs, yet operate on different data configurations. Presence-only models only require the present evidence of the modelled phenomenon, while the counterpart group, presence-absence models, would require both presence and absence locations to better distinguish the presence patterns from absences. When it comes to modelling natural phenomena, selecting the absence locations may become troublesome. It is because some areas may look deprived of any morphological appearance of the phenomenon, yet they may be subjected to different degrees of susceptibility/potentiality. In other words, it is difficult to categorize an area as an absolute absence location. Hence, presence-only models gain the upper hand in this specific matter.
Maximum entropy is a widely used presence-only model that inherits a robust concept based on which maximizing the entropy helps the model better distinguish the  Phillips et al. (2004Phillips et al. ( , 2006, a versatile platform is bestowed upon the end-users by many beneficial options. In a nutshell, MaxEnt starts extracting the values of all controlling factors at the areas that the evidence is present (stored as point features). Then, it transforms the mathematical form of the factors into another set of arithmetic functions, called features (Jaafarzadeh et al. 2021; Mousazade et al. 2019). Further, the model takes a first arbitrary guess that the probability of observing the phenomenon (or the set of conditions that can lead to the presence of the phenomenon) in all pixels equals 0.5. As the process of maximizing the entropy using the Gibbs distribution function and Lagrangian solutions continues, the first guess becomes more and more precise, and the first probability distribution function (pdf) becomes more similar to the target pdf (i.e., actual pdf of the phenomenon of interest). As mentioned above, two sets of data are required to run the MaxEnt model: evidence, which is the location of areas that represent the high potential of groundwater (wells and springs), and controlling factors (continuous and categorical types such as slope and land use/cover). For the latter, we strived to extract the most relevant morphometric indices from the DEM layer, including slope, aspect, convergence, flow accumulation, negative and positive openness (surface convexities and concavities from a given zenith), plan and profile curvature, relative slope position (RSP), stream power index (SPI), topographic position index (TPI), terrain ruggedness index (TRI), and topographic wetness index (TWI). All these factors are considered topo-hydrologic drivers that can govern groundwater potential. Together with the interpolated double-ring tests conducted during the field surveys, the above-listed factors established the continuous data. As for the categorical data, the land use/cover map was acquired from the organizational sources and updated based on Google Earth images and field visits. Thematic maps used for groundwater potential modelling are presented in Fig. 4.
Regarding the evidence layer, the well locations were the only representative data available of groundwater potential. The insufficient number of wells would inevitably affect the reliability of the data mining model due to the large extent of the study area and the necessity of data mining models to provide them with big data. Therefore, a well-established data enrichment technique was taken into account. Drawing on a study conducted by Tsangaratos and Benardos (2014), Kornejady et al. (2017), and Davoudi Moghaddam et al. (2020), the MD was found to be a valuable technique for data enrichment and, accordingly, boosting the performance of machine learning models. The MD technique operates on the inverse covariance matrix of the controlling factors and considers their spatial dependencies, following Eq. 5 (Fukunaga 1990).
where MD denotes the Mahalanobis distance, x is the vector of each factor, m is the mean value of each factor, T is the transposed matrix of the factors' values, and ∑ −1 is the inverse covariance matrix.
We applied the MD technique in ArcGIS 10.8 using the continuous and categorical data listed previously and the 21 well locations. The derived MD map with distance values, in which the smaller distances represent conditions similar to high groundwater potential, was rescaled on a 0-1 scale using the chi-square method (the number of conditioning factors minus 1 was treated as the degree of freedom). The final map was then classified into two groups using the natural break classifier in which the group with higher values within (close to 1) was selected to enrich the presence data. A total of 1000 locations were drawn and treated as presence locations (i.e., evidence layer) for the MaxEnt model.
A boosting technique with 500 iterations, a convergence value of 10 −5 , 10,000 background points (i.e., pseudo absence locations) and auto-feature (i.e., automatic feature selection) were used as the standard configuration of Max-Ent and a rule of thumb proposed by Phillips et al. (2004Phillips et al. ( , 2006. Furthermore, due to the 1000 enriched presence data, a 50:50% data sampling method was adopted to split the data into two sets of training and validation (randomly partitioned in MaxEnt software), guaranteeing data noise removal. It is noteworthy that the sample partitioning method is not fully random but rather adopts a systematic-random approach. In a nutshell, the MD map can be rendered on a 0-1 scale, where values close to 0 represent higher distance values; that is, higher distances from groundwater samples logically imply 0 or close to 0 groundwater potentiality. In addition, we split the 0-1 scaled map into two portions: 0-0.5 and 0.5-1. As previously mentioned, the latter is used to randomly generate 1000 samples because they are representative of moderate to high groundwater potentiality. Now we can randomly split these 1000 points into 50% partitions and use one half for training and another for validation. However, in order to train the MaxEnt model most suitably, we adopted a systematic-random approach. In it, once again, we split the 0.5-1 into two portions: 0.5-0.75 and 0.75-1, and randomly selected the 50% training dataset in such a way that 25% of the samples fall within 0.5-0.75 and the other 25% in 0.75-1 portion. The same strategy was applied to select the 50% validation dataset. In this way, the model trains based on the sample points that are spatially well scattered. Similarly, the generalization capacity of MaxEnt is validated in a more representative manner. If this strategy was not taken into account, the possibility of randomly selected samples concentrated in a small area will most likely remain to exist, which can deteriorate the generalization/prediction power of MaxEnt and also the credibility of the entire modeling process (Rahmati et al. 2016). MaxEnt model performance was assessed using the receiver operating characteristic (ROC) curve, plotted using both training and validation datasets. The area under the ROC curve (AUROC) is a cutoff-independent metric that reflects the goodness-of-fit and prediction power of the model once, respectively, fed by training and validation samples. The x-axis plots the false-positive error (i.e., locations incorrectly predicted as high groundwater potential by the model) against the y-axis as the true positive (i.e., locations correctly predicted as high groundwater potential by the model, in line with field observations). The AUROC values close to 1 represent a well-performing model, while values close to 0.5 signify a neutral model, completely acting out of random chance (Pontius and Schneider 2001).
The following equations were used to complete the ROC curve:  where TP represents the true positive, FN represents the false negative, TN represents thee true negative, FP represents the false positive, P indicates the number of total gullies, and N indicates the number of total non-gullies (Arabameri et al. 2020). Table 2 details the areal extent of the vulnerability classes based on the DRASTIC, GODS, and SI models. Based on the DRASTIC model, the vulnerability map of the aquifer was divided into three categories of low (47-92), moderate (93-136), and high vulnerability (137-184). Based on this classification, the moderate and low vulnerability classes represented the highest and lowest portion of the study area. Based on the results of the GODS model, vulnerability zonation of the aquifer was classified into two classes of low (0.1-0.3) and moderate (0.3-0.5), in which the low vulnerability class accounts for the largest area in the region. Moreover, the SI method divided the aquifers into two classes with almost equal extents. Figure 5 depicts the final vulnerability maps derived from DRASTIC, GODS, and SI.

The validity of the vulnerability maps derived from the SI, GODS, and DRASTI models
Nitrate concentration data of 21 wells in the region were used to validate the SI, GODS, and DRASTIC models. Figure 6 shows the location of 21 well used to determine nitrate concentration and the validation process. By superimposing the nitrate concentration layer on the three adopted models' final vulnerability maps, all the points with average nitrate concentration fall within the district with moderate contamination vulnerability of GODS and SI models. That is, the districts with average nitrate concentration accurately exhibit moderate vulnerability. This result indicates that the accuracies of these models are satisfactory. The districts with low nitrate concentration also match with the districts with low nitrate vulnerability class. In general, a comparison of nitrate concentration and vulnerability maps of different models showed that the DRASTIC model had a higher consistency and accuracy than SI. Besides, the GODS model's accuracy was less than the two other methods such that in some districts of the aquifer, highly vulnerable regions reported by the DRASTIC model had high nitrate concentration.
According to the field experiments of the NO −3 level at the well locations and plotting the vulnerability values derived from SI, GODS, and DRASTIC models against the nitrate values (Table 3), a simple yet practical analysis of the validity of groundwater contamination models was feasible. The linear regression equations and the corresponding R 2 values are presented in Fig. 7. Based on the R 2 values, it is evident that the vulnerability values derived from the DRASTIC model are strongly correlated with nitrate concentrations, making the DRASTIC the superior model for assessing the vulnerability of the aquifer to contamination sources in the Bahabad region.

Groundwater potential mapping using the coupled MaxEnt-MD model
A total of 1000 training and validation sample points were extracted from the rescaled MD model (Fig. 8) and fed to the MaxEnt model along with the 13 thematic layers of the controlling factors. The final product of the MaxEnt model was then classified into three potential groups using the natural break classification scheme (Fig. 9). The corresponding ROC curves of the training and validations stages of the MaxEnt model are presented in Fig. 10, indicating a magnificent performance of MaxEnt in groundwater potential modelling. Blue, green, and red bars in the jack-knife test presented in Fig. 11, respectively, indicate the AUC of the model using each factor alone (i.e., considering each factor individually), one-by-one factor removal (i.e., considering all the factors but one), and using all the factors. Higher AUC values in the blue bar represent the high contribution of each specific factor to the modeling process. The results revealed that elevation, positive openness, and RSP contribute the most to groundwater potentiality in the study area, followed by slope, TRI, negative openness, and infiltration rate (doublering test). On the other hand, the lower AUC in the green bar represents the missing unique information due to the removal of a specific factor. In this regard, elevation and infiltration rate exhibit the most unique information for Max-Ent. Based on the overall AUC (red bar), it is evident that using all the factors would gain the highest AUC value and therefore a better-trained model.

Discussion
In the DRASTIC method, moderate vulnerability occupied the largest area in the region, while high vulnerability occupied the smallest area. A part of the northern parts of the aquifer (0.05%) is in the area with high vulnerability, most of the region (91.23%) is in the area with moderate vulnerability, and some scattered areas in the eastern, western, and southern parts of the region (8.72%) face low vulnerability. According to the GODS model, the southern and northern parts of the plain (ca. 60%) have a low vulnerability, while areas with medium vulnerability (ca. 40%) are located in the central parts. According to the SI method, in terms of vulnerability, about 50.45% of the region is in the low vulnerability range, while the rest of the region is in the moderate vulnerability range. Based on the results obtained from the three models and their corresponding maps, the central and southern parts of the plain are more vulnerable than other areas, which include pasture and agricultural uses. In general, the DRASTIC model provided higher accuracy than the other two models.
Based on the results of the vulnerability map derived from the DRASTIC model, most areas in the east, west, and south had moderate vulnerability, while high vulnerability is concentrated in small portions of the region. The GODS method showed that eastern, central, and Estern and some small parts of the north exhibit moderate vulnerability, while northern and southern parts of the region have low vulnerability. The SI model gives a somewhat more discriminative vulnerability map where low and moderate vulnerability regions subsequently appear contiguously. The differences in spatial patterns of low-to-high aquifer vulnerability between the models stem from their different algorithms for identifying the main hotspots.
The main benchmark to analyze the performance of vulnerability models to contamination is the examined nitrate concentration. Accordingly, the DRASTIC model is discernibly well-performed and has a statistically significant and strong positive correlation with nitrate concentration. The main reason behind the success of the DRASTIC model can only stem from being more inclusive and incorporating more groundwater quality-representative and -relevant factors into its computational algorithm. Literature review attests to these inferences, where DRASTIC has been previously found to be a better tool to assess groundwater vulnerability to pollutants due to conglomerating multiple data from the hydrosphere, lithosphere, and anthroposphere even for multiple aquifers (Nadiri et al. 2017). Moreover, Zhao and Pei (2012) investigated groundwater contamination using different models, including GODS, DRASTIC, and GLA, in a study site in China. They reported that among the adopted models, the DRASTIC model performed better than the others, and the inclusion of more parameters has been reported as the main plausible reason for this success. Also, Ahmadifar et al. (2017) investigated the contamination of groundwater resources using DRAS-TIC and SINTACS in Sarab plain, Iran, and reported that the DRASTIC showed better results than the SINTACS model. Stefania et al. (2019) carried out groundwater contamination assessment in Italy using various algorithms and mathematical models. The results of their research also emphasize the vital role of incorporating more relevant parameters in gaining better and more reliable results.
A simple visual check of the derived vulnerability maps from the corresponding models indicates spatial asymmetries between their results. This notion can stem from the different computational algorithms of the models; however, it can be alarming for end-users, and a more robust and rigorous framework for all the modelling stages, be it model selection, the configuration of inputs, or performance assessment, would be critical. The strong correlation between the nitrate level and the vulnerability values obtained from the DRASTIC model attests to its applicability to modelling the groundwater vulnerability to contamination sources.
Overall, the inferences offered by groundwater contamination models are required for making reliable decisions on proactive measures for complex aquifer systems (Nadiri et al. 2017). In particular, in arid areas with low water inflow and excessive water consumption from aquifers, and the presence of polluting sources can make it difficult for the aquifer systems to withstand the pressing meteorological and anthropogenic agents and stresses, which adds to the importance of modelling groundwater vulnerability to contaminants. This matter should be taken into account existing opportunities and suppressing threats. The latter was carried out in this study by modelling the groundwater potential mapping.
The ROC curves plotted for the MaxEnt model attested to its magnificent performance in extracting the latent pattern of groundwater potential in the region, which is in line with previous studies (Rahmati et al. 2016;Razandi et al. 2017;Raney and Leopold 2018). The high value of training performance (i.e., 0.907) indicates a high goodness-of-fit and an almost perfect learning stage thanks to 500 iteration sessions. Moreover, the overfitting issue was obviated by the high AUROC value for the validation stage (that is, 0.901), which attests to the high prediction power and generalization capacity of MaxEnt.
Due to the limitation of available data and the lack of distribution of data for the entire domain level, we had to use the Mahalabois distance (data enrichment). In the Mahalanobis interval, outlier data are also detected and normalized, and the correlation between factors is considered.
The jack-knife test revealed the high contribution of elevation to the modelling of groundwater potential. The saturation runoff generation mechanism (wetlands/ waterlogged) is mainly seen in low altitudes, where the groundwater table approaches the Earth's surface. Higher values of positive openness indicate convexity, where water can store and gradually percolate into the deep strata, which gives a connotation of a high groundwater potential. The RSP (relative slope position) also delineates valleys and flat areas with high groundwater potential from high ridges, thus contributing greatly to the modelling process.
The importance of elevation as a proxy factor to multilateral topo-hydrological attributes was evident in the sensitivity analysis in the jack-knife test. The removal of infiltration rate led to a decrease in the AUC value of MaxEnt, indicating its important role in identifying the high groundwater potentiality. Areas with a high infiltration rate are mostly The groundwater potential map obtained from the MaxEnt model and the corresponding potential classes covered by coarse-grained materials with high porosity that would expedite permeability.
The interactive analysis between the vulnerability map derived from the DRASTIC model and the groundwater potential map derived from MaxEnt went beyond and scrutinized the spatial impact of contamination sources. In this sense, the poor correlation can somewhat obviate local concern regarding groundwater contamination and implies that although a large part of the area is exposed to moderate vulnerability, areas with high groundwater potential are still safe and are located far from the high vulnerability zone. However, this can be a blueprint for the authorities to take the contamination sources more seriously and lean towards interactive and multidisciplinary approaches in order to gain a holistic view and a trade-off between the prevailing threats and opportunities upon the region.

Conclusion
Accurate modelling and zoning of groundwater contamination vulnerability in arid areas with limited water resources are of great importance. The general inferences of the adopted experimental models (i.e., DRASTIC, GODS, SI) in the Bahabad region shows that a different combination of aquifers' physical attributes can lead to an entirely different pattern of groundwater vulnerability. The DRASTIC model, compared to its counterparts, incorporates more Fig. 10 The ROC curve plotted for the training and validation (test) stages

Fig. 11
Jack-knife test obtained from MaxEnt model, indicating the AUC values in three forms of solitary factor adoption (blue bar), one-by-one factor removal (green bar), and engaging all the factors (red bar) (aspect: slope aspect, dem: elevation, doublering: infiltration rate, flow_acc: flow accumulation, lulc: land use/cover, no: negative openness, po: positive openness, rsp: relative slope position, slope: slope degree, spi: stream power index, tpi: topographic position index, tri: terrain ruggedness index, twi: topographic wetness index) data and gives a more representative result, which was confirmed by the nitrate concentrations examined in 21 wells. As appears, most of the area has moderate vulnerability to contaminants. Considering the transient status of Iran's plain and imminent natural disasters such as droughts, this should be taken seriously and act as an alarming status. However, the groundwater potential pattern modelled by MaxEnt shows that there is still a chance to recover. More precisely, the current groundwater potential does not coincide with the high vulnerability status. The application of MD ensemble with a machine learning model showed promising results for data enrichment, which should be considered in areas with scarce data. This work promotes the reciprocal groundwater potential and contamination modelling (GPCM), which is currently neglected by the trail of studies disseminated by the scientific community and the management sectors in Iran.
Our upcoming research conducts a more detailed analysis of the geochemical tests conducted at the operating wells in the study area and supports the high performance of the DRASTIC model in a more elaborative manner.
The slope factor plays an important role in determining the velocity with which the water flows. Slope controls the infiltration of groundwater into subsurface and recharge processes (Adiat et al. 2009). Slope and influence have opposite relationship. In such a way that on a gentle slope, the surface runoff moves slowly and has more time to penetrate, and steeper slopes have more runoff and the underground water penetrates less. The direction of the slope plays a role in the amount of rainfall, the effects of wind and receiving sunlight. As the altitude increases, the groundwater potential decreases. The reason for that is the formation of aquifers with a large thickness at high altitudes. Also, in high altitudes, due to the increase in slope, more runoff is produced and the amount of nutrition is reduced. By increasing the value of the topographic position index (TPI), it indicates high and hilly areas, which reduces the underground water potential. An increase in SPI causes an increase in mud pollution and, as a result, a decrease in infiltration (Golkarian and Rahmati 2018). By increasing the position of the relative slope (RSP), we have a decrease in the underground water potential. The use of pastures, gardens, residential areas, and agriculture due to excessive consumption has a great impact on the underground water potential and reduces it (Teimouri and Asadi Nalivan 2021).
One of the limitations of this project is the huge cost of information and the lack of input information, and it is necessary to mention that the quantitative and qualitative information of the area was very little and limited, but due to the importance of the study area, this area was used as research area.
It is suggested that quantitative models such as mudflow and white box models should be used in future studies by spending a lot of money and collecting well points, so that these models can check uncertainty and error checked more accurately than qualitative models.