Mapping of Groundwater Spring Potential in Karst Aquifer System Using Novel Ensemble Bivariate and Multivariate Models

Abstract

Groundwater is an important natural resource in arid and semi-arid environments, where discharge from karst springs is utilized as the principal water supply for human use. The occurrence of karst springs over large areas is often poorly documented, and interpolation strategies are often utilized to map the distribution and discharge potential of springs. This study develops a novel method to delineate karst spring zones on the basis of various hydrogeological factors. A case study of the Bojnourd Region, Iran, where spring discharge measurements are available for 359 sites, is used to demonstrate application of the new approach. Spatial mapping is achieved using ensemble modelling, which is based on certainty factors (CF) and logistic regression (LR). Maps of the CF and LR components of groundwater potential were generated individually, and then, combined to prepare an ensemble map of the study area. The accuracy (A) of the ensemble map was then assessed using area under the receiver operating characteristic curve. Results of this analysis show that LR (A = 78%) outperformed CF (A = 67%) in terms of the comparison between model predictions and known occurrences of karst springs (i.e., calibration data). However, combining the CF and LR results through ensemble modelling produced superior accuracy (A = 85%) in terms of spring potential mapping. By combining CF and LR statistical models through ensemble modelling, weaknesses in CF and LR methods are offset, and therefore, we recommend this ensemble approach for similar karst mapping projects. The methodology developed here offers an efficient method for assessing spring discharge and karst spring potentials over regional scales.

1. Introduction

Groundwater has been considered the most important resource in arid and semi-arid areas. Population increases and economic growth have resulted in water shortages around the globe, especially in developing countries, leading to greater exploitation and reliance on groundwater resources [1,2]. Among the various forms of subsurface freshwater, karst groundwater is a major resource for water supply in many regions, including China, Turkey, Iran, Mexico, and the United States [3,4,5,6,7]. Karst aquifers are relatively high-yielding, heterogeneous systems that are characterized by complex groundwater flow paths [8,9,10]. They store and supply drinking water for approximately 20–25% of the world’s population, while constituting only small areas of the global land surface [11]. The characterization of karst aquifers is hampered by their complex character, and the often limited number of wells available to support hydrogeological observations and geological descriptions [12].

Groundwater models of karst aquifer systems invariably simplify karst features—and the spatial heterogeneity more generally—because it is impossible to acquire sufficient data to properly delineate the main flow features [13]. Additionally, specialized models are required to simulate the conduit flow that occurs within karst aquifers, at least for situations where equivalent porous medium approaches are not valid [14]. However, physically based models of karst aquifers are limited by both data availability and the heavy computational burden of conduit-matrix simulation. Thus, strategies to characterize the hydrogeology of karst aquifers, and in particular, to quantify spring hydrology, are essential precursors to the effective management of groundwater in these systems [15].

Geographical information systems (GIS) and remote sensing (RS) are commonly used to evaluate the hydrogeology of karst aquifers, including the spatial distribution of spring potential [16,17]. Recently, Chen et al. [18] undertook a global-scale project—the World Karst Aquifer Mapping (WOKAM)—to map karst aquifers and carbonate rock settings around the world. Karst springs, wells, caves and other water abstraction structures were also identified, although the large scale of the project meant that many regions have only limited information regarding point features within the landscape.

‘Spring groundwater potential’ was defined as the possibility of spring occurrence (with sufficient water) in an area and helps planners to estimate groundwater productivity [19,20]. There are myriad approaches to groundwater potential mapping using GIS-based statistical models. These include: frequency ratio (FR) [16,21], weight-of-evidence [22], evidential belief function (EBF) [23], logistic regression (LR) [24,25], certainty factor (CF) [26], analytical hierarchy process (AHP) [27], Shannon’s entropy [28], maximum entropy [29], support vector machine (SVM) [30] and boosted regression tree [31]. Groundwater potential mapping highlights the locations of potential groundwater zones including low, moderate and high for the optimal use of resources [32]. In general, multivariate statistical models (e.g., LR) cannot analyze the classes of independent factors (e.g., lithology, land use/land cover, soil, etc.) because they just determine the weight of predictive factors. In addition, bivariate statistical models (e.g., CF) are limited in their capacity to recognize the importance of independent variables [26] because they use a class-based procedure in modelling. To overcome these problems, scientists have developed ensemble statistical approaches [33,34,35,36,37,38,39].

Ensemble modelling is the combination of multiple predictive models that are used conjunctively to provide an evaluation of the same dataset [39]. Ensemble modelling approaches have been widely used in different fields, including groundwater [40,41,42,43,44,45,46,47], flood [40,48,49,50,51,52,53,54,55,56,57], landslide hazard [58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87], land/ground subsidence [88,89], gully erosion [90,91], dust storm [92], wildfire [93], sinkhole [94], droughtiness [95,96], earthquake [97,98] and species distribution [99,100,101]. These studies provide clear evidence that the application of ensemble models can potentially result in improved capability (over individual methods) of GIS-based statistical models [102]. However, previous studies aimed at the prediction of karst spring occurrence have not investigated the potential of ensemble models that integrate bivariate and multivariate statistical models to support the construction of spring potential maps [24,25].

The main aim of this study is to explore different statistical models, both single and ensemble, to assess spring groundwater potential in a karst area, namely Bojnourd Region, northeast Iran. Although various researchers have tried to delineate groundwater spring potential zones [30,45,103,104], the evaluation of other modeling approaches is highly necessary. To address this, a new ensembling technique is proposed to improve the accuracy of the statistical models and to identify groundwater potential zones. The novel approach of ensemble certainty factor and logistic regression models has not previously been touched upon in the literature. Therefore, the principal motivation of this study is to use this ensemble model to delineate prone areas to karst springs in Bojnourd region, Khorasan Shomali, Iran. It is anticipated that the outputs of this study will identify the strengths and weaknesses of the applied models, thereby, providing improved insights into model choice in the spatio-temporal analysis of karst springs discharge using statistical models.

2. Study Area

The study area, Bojnourd Region, is located in the north east of Khorasan Shomali province, Iran, between latitudes 37°15′ to 37°35′ N and longitudes 57°03′ to 57°40′ E (Figure 1), with minimum and maximum elevation of 875 and 2968 m, respectively. The mean precipitation is 272 mm/year and temperature ranges between 6.8 to 19.7 °C [105]. There are six main lithological formations in the study area including Aitamir (Cretaceous), Sarcheshmeh (Early Cretaceous), Shurijeh (Jurassic-Cretaceous), Upper-Red (Miocene), Niur (Silurian) and Shemshak (Triassic-Jurassic). Lithological characteristics of these formations are summarized in

Table 1. Geological formations in the study area.

Geology (Group)	Formation	Lithology
Cretaceous (C)	Aitamir	Olive green glauconitic sandstone and shale
Early Cretaceous (EC)	Sarcheshmeh	Ammonite bearing shale with interaction of orbitolin limestone
Jurassic-Cretaceous (JC)	Shurijeh	Pale red argillaceous limestone, sandstone and conglomerate
Miocene (M)	Upper-Red	Red marl, gypsiferous marl, sandstone and conglomerate
Quaternary (Q)	-	Low level piedmont fan and valley terrace deposits
Silurian (S)	Niur	Coral limestone and dolomite, shale, sandstone
Triassic-Jurassic (TJ)	Shemshak	Subordinate sandy limestone, dark grey shale and sandstone

. From a geomorphological point of view, the study area is dominated by a mountainous landscape. In addition, there are numerous folds and faults in the study area which cause the occurrence of springs [106]. The major faults are oriented NE-SW and the others NW-SE and N-S. Inceptisols and entisols are dominant soil types of the study area. The depth of soil varies from one to three meters, although there are numerous surficial rocks (lithological outcrops) in steep areas of upstream zones. According to the reports of Iranian Department of Water Resources Management (IDWRM), hydraulic conductivities in this mountainous region range between 10⁻⁴ and 10⁻⁵ m/s.

Nearly 88% of water demand in Khorasan Shomali province is supplied by groundwater resources [107]. Therefore, it has become necessary to predict groundwater potential zone in this region. The discharge of springs ranges between 0.5 and 15.2 SI. The size of karstic areas in this province is over 17,000 km², which are used as karst water resources to provide drinking water for the Khorasan Shomali and Khorasan Razavi provinces. Furthermore, Khorasan Shomali province covers approximately an area of 248,000 km², which is inhabited by more than 6 million people and where agriculture is the major activity.

3. Methodology

The adopted methodology (Figure 2) included: preparation of data, construction of models and production of karst spring potential map (KSPM), validation of results and evaluation of their efficiency.

3.1. Data Used

In order to analyze karst spring potential, springs locations and 10 spring-affecting factors were considered as dependent and independent variables, respectively. All data were obtained from the Iranian Department of Water Resources Management (IDWRM).

3.1.1. Karst Spring Inventory

A total of 359 karst springs was identified and recorded by the Iranian Department of Water Resources Management (IDWRM) in the study area, of which 251 (70%) were randomly selected for calibration/training of models and the remaining 108 (30%) for validation purpose. The distribution of training and validation karst springs is shown in Figure 1.

3.1.2. Karst Spring-Affecting Factors

As mentioned by Crozier [108], preparation of corresponding thematic data layers and the selection of factors are important components of any model. In this study, 10 karst spring-affecting factors including distance from fault, land use/land cover, geology, topographic elevation, slope angle, slope aspect, plan curvature, topographic wetness index (TWI), drainage density and distance from river were utilized for producing KSPMs.

Distance from fault

Distance from fault is among those factors that has been widely used as an important groundwater influential factor [109]. The faults were extracted from a geological map of the study area as line features in ArcGIS 10.2 and a map of ‘distance from faults’ was developed using the Euclidean distance tool and classified (Figure 3a).

Land use/ land cover

Land use/land cover is a significant factor in controlling groundwater recharge processes [110]. The land use/land cover map of the study area (Figure 3b) shows that rangeland covers most of study area, followed by agriculture, forest and urban classes (Figure 3b).

Geology

Since the hydraulic conductivity and porosity of soils and rocks are affected by lithological variations, groundwater occurrence is highly influenced by geology [104,111]. In this study, the geologic map of study area with seven lithological groups was extracted from 1:100,000 scale geologic map of Iran (Figure 3c, Table 1).

Topographic factors

Surface topography controls the runoff rate and direction [112,113]. Therefore, topographic elevation has been frequently used by researchers as a controlling parameter for groundwater occurrence [31,114]. The digital elevation model (DEM) map of the study area was obtained from the IDWRM and we classified it using an equal-interval classification scheme (Figure 3d). IDWRM originally produced the DEM from aerial photographs and topographical lines was obtained from the Iranian Department of Water Resources Management (IDWRM) at 20 m cell resolution. Slope angle affects the infiltration, subsurface flow and duration of overland flow [115,116]. Therefore, the groundwater recharge process is strongly controlled by slope angle [117]. According to Magesh et al. [115], due to low infiltration and high run-off in steeper slope, these areas have poor potential to groundwater accumulation. The slope angle map of study area was extracted from the DEM and classified into seven classes based on IDWRM classification (Figure 3e). Plan curvature represents the rate of slope aspect or slope gradient change. Negative and positive values of plan curvatures characterize concavity and convexity in the downslope direction, respectively. Zero curvature also represents flat areas [118]. As mentioned by Lee and Pradhan [118], concave slope probably contains more water. The plan curvature map of study area was produced from DEM in concave, flat and convex classes (Figure 3f). The slope aspect layer was derived from DEM and divided into nine classes including eight main directions (Figure 3g).

Water-related factors

There is a significant relationship between distance from river and groundwater spring occurrence [31]. The maps of drainage density and distance from rivers were created from the topographic data using Kernel density and Euclidean distance tools, respectively. The maps of drainage density and distance from river are shown in Figure 3h,i, respectively. The soil moisture conditions are reflected in values of TWI [119], which is defined as:

$$TWI=ln\left(\frac{A_s}{tan\beta }\right)$$

(1)

where A_s is the specific catchment area (m²/m) and β is slope gradient [120]. TWI map was shown in Figure 3j. Drainage density is defined as ratio of total length of a stream to the area of the drainage basin [121]. In general, low drainage density areas are more suitable for groundwater development in comparison with high drainage density [122].

3.2. Statistical Methods

A bivariate statistical analysis (certainty factor) and a multivariate statistical analysis (logistic regression) were applied individually and as an ensemble to delineate karst spring potential. Finally, the ensemble map was compared with each individual map using the receiver operating characteristic (ROC) analysis to find if there is any improvement in combining the models.

3.2.1. Certainty Factor Model

Shortliffe and Buchanan [123] proposed and defined certainty factor (CF) as a probability function that was later modified by Heckerman [124]. Different researchers have applied this model in their studies for landslide susceptibility mapping [125] and groundwater potential mapping [26]. The formula to calculate CF is [124]:

$$CF=\{\begin{array}{c}\frac{P_a-P_b}{P_a\left(1-P_b\right)} if P_a\ge P_b\\ \frac{P_a-P_b}{P_b\left(1-P_a\right)} if P_a<P_b\end{array}$$

(2)

where P_a is conditional probability of having a number of karst spring occurrence in class a, P_b is prior probability of having the total number of karst spring occurrence in the study area. The CF value for different attribute classes in each conditioning factor were calculated from their relationship to the spring occurrence. The CF values ranges from −1 to +1. A positive value of CF means that the groundwater potential is higher, while a negative value illustrates a decrease in the certainty of groundwater occurrence [126]. In order to delineate the combination of two different layers of information, which is denoted as Z, the CF values of the spring-affecting factors are pairwise combined using the integration rules as follow [127].

$$Z=\{\begin{array}{cc}C{F}_X+C{F}_Y-C{F}_{X}C{F}_Y& C{F}_X and C{F}_Y\ge 0\\ \frac{C{F}_X+C{F}_Y}{1-min\left(\left|C{F}_X\right| and \left|C{F}_Y\right|\right)}& C{F}_{X}and C{F}_Y Opposite sign\\ C{F}_X+C{F}_Y+C{F}_{X}C{F}_Y& C{F}_X and C{F}_Y<0\end{array}$$

(3)

where CF_X and CF_Y are CF values of two different spring-affecting factors.

The pairwise combination is carried out in different iterations until all the spring-affecting layers are added to prepare the KSPM. The final map was categorized into four classes: low, medium, high and very high [29,128].

3.2.2. Logistic Regression Model

LR aims at finding the best-fit model to express the relationship among a binary dependent variable (groundwater occurrence) and some independent variables [129]. This model is helpful to predict the absence or presence of outcome based on predictive variables (Lee 2005). The LR model is widely applied in different fields of geosciences [24,130]. However, it has some limitations to determine the weight of classes of independent variables [131]. Hosmer and Lemeshow [132] and Kleinbaum et al. [133] also explained the detailed information of this technique. The general form of logistic regression is as follows:

$$Z=b_0+b_1{x}_1+b_2{x}_2+\dots +b_m{x}_m$$

(4)

where Z is a linear combination function of the explanatory variables, x₁, x₂, …, x_m are independent variables and b₁, b₂, …, b_m are the coefficients of the logistic regression to be estimated [134]. The relationship between the occurrence and its dependency on several variables is defined as:

$$P=\frac{e^z}{1+e^z}$$

(5)

where P and e are the probability of spring occurrence (varying from 0 to 1) and the Neperian number, respectively [135]. In addition, Z factor was defined above Equation (4). Value of 1 means the presence of a spring, and value 0 indicates the absence of a spring [24]. After creating the KSPM by LR model, we classified it using natural break classification method. As mentioned by Hosmer and Lemeshow [132], independent variables in LR model is sensitive to co-linearity for model fitting, while more than two variables in LR model is often called ‘multi-collinearity’. Multi-collinearity refers to a linear relation among two or more variables [136]. LR model was used to do statistical analysis and determine whether there was a significant effect on spring occurrence. Tolerance (TOL) and the variance inflation factor (VIF) are significant indicators for determining multi-collinearity among the independent variables [137] that were applied (

Table 2. Indices of multi-collinearity diagnosis for independent variables used in logistic regression (LR) modelling.

Spring-Affecting Factors	TOL	VIF
Plan curvature	0.89	1.18
Slope aspect	0.96	1.20
Slope angle	0.86	1.05
Distance from river	0.94	1.20
Drainage density	0.84	1.08
TWI	0.83	1.17
Topographic elevation	0.93	1.06
Geology	0.97	1.13
Land use/ land cover	0.87	1.02
Distance from fault	0.82	1.05

). Based on this diagnosis method, variables with TOL < 0.1 and/or VIF > 10 shall be excluded from further analysis [138]. Moreover, the overall statistic of the LR model for all 10 independent variables was assessed using χ² value of the Hosmer–Lemeshow test, Cox and Snell R², Nagelkerke R² and pseudo-R² mesaures.3.2.3. Ensemble model.

Ensemble modeling is a technical process of synthesizing the results of single models into a single incorporated model for enhancing the prediction [139]. This type of modeling approach tends to overcome their individual disadvantages while keeping their individual advantages [102]. Weighted linear combination technique is proposed here as an ensemble method to combine the result of CF and LR:

$$\begin{array}{c}KSPM=\left(D{F}_{r}D{F}_w\right)+\left(L_r{L}_w\right)+\left(G_r{G}_w\right)+\left(A_r{A}_w\right)+\left(TW{I}_{r}TW{I}_w\right)+\left(D{D}_{r}D{D}_w\right)\\ +\left(D{R}_{r}D{R}_w\right)+\left(S_r{S}_w\right)+\left(A{S}_{r}A{S}_w\right)+\left(P{C}_{r}P{C}_w\right)\end{array}$$

(6)

where DF is distance from fault, L is land use/ land cover, G is geology, A is topographic elevation, TWI is topographic wetness index, DD is drainage density, DR is distance from river, S is slope angle, AS is slope aspect, PC is plan curvature and the subscripts ‘r’ and ‘w’ indicate ratings and weights for spring-affecting factors, respectively. Ratings are obtained from the CF model (i.e., the weights of classes of a given spring-affecting factor), whereas w was calculated based on the spatial sensitivity analysis (SSA) of LR model. The SSA allows to determine the influence of predictive variables on the model output [16,140]. In this study, this was conducted by excluding each spring-affecting factor and determining w through the relative decrease of AUC values Equation (7) [29]:

$$w=\frac{\left[AU{C}_{all}-AU{C}_i\right]}{AU{C}_{all}}\times 100$$

(7)

where AUC_all and AUC_i indicate the AUC values obtained from the LR model using all spring-affecting factors and the prediction when the ith spring-affecting factor was excluded, respectively.

3.2.3. Performance Assessment

Performance assessment is required in any predictive analysis [141]. Nampak et al. [128] also stated that there is no merit in modeling with no scientific significance; hence, validity and uncertainty of models need to be assessed. Furthermore, many researchers have validated their results using different models. Naghibi et al. [31], Pourghasemi and Beheshtirad [109] and Rahmati et al. [3] have utilized area under ROC curve (AUC) to assess the accuracy of groundwater potential maps. This method illustrates to what extent the occurrence or non-occurrence of a considered occurrence is correctly predicted [29]. In the ROC curve, false and true positive rates are plotted in X- and Y-axis, respectively. The resultant AUC value ranges from 0.5 (i.e., a random prediction) to 1 (i.e., an excellent prediction) [142].

4. Results

4.1. CF Modelling

CF values of spring-affecting factors are presented in Figure 4. As depicted in the figure, negative and positive relationships between factors and spring occurrence are well illustrated. Obviously, negative relationship exists between spring potential and distance from faults parameter. Areas with more distance from fault have lower potential to spring occurrence. In the case of geology factor, some lithologic formations including Silurian (S), Jurassic-Cretaceous (JC) and Triassic-Jurassic (TJ) had positive influence in spring occurrence with CF values of 0.94, 0.42 and 0.21, respectively.

In the case of land use/land cover, urban areas and agriculture have the least and most impact on groundwater, respectively. The impervious surface areas in urban areas tend to reduce recharge and hence groundwater availability. As the agriculture areas as a consequence factor need to high water requirement, these areas have a higher potential to groundwater recharge [47]. It was also found that there was no significant correlation between spring occurrence and topographic elevation classes. Areas with elevation ranges of 1676–1978 m and 1414–1676 m have positive influence on spring occurrence probability with CF values of 0.45 and 0.16 which indicate a good karst spring occurrence probability. In another study in this region, Rahmati et al. [44] indicated that spring occurrence probability will reduce as elevation decreases which is in agreement with the results of our study. It is noted that, groundwater is more probable to occur in areas with high TWI values because it is well apparent in Figure 4e that positive correlation exists between increasing spring occurrence and high TWI value.

Conversely, Razandi et al. [26] stated that there is not positive correlation between groundwater occurrence and TWI, which shows a higher groundwater potential over a decreasing TWI value. The reason of this difference between their results and our achievements refer to geomorphological characteristics of study areas (i.e., their study area were plain whereas our study area is a mountainous region). In areas where drainage density ranges from 1.3 to 2.8 km/km² (CF = 1.82), groundwater potential is higher indicating that springs are more probable in areas with denser drainage density. The assessment of distance from river demonstrated that areas with lower distance from river have more potential for spring occurrence and as this distance increase, the probability to groundwater occurrence decreases. Lee et al. [143] proved these results by indicating that drainage density and distance from the river have positive and negative influence on groundwater yield. Results indicated that the areas with 8–40° slope angles have positive influences with highest value in 14–20° and 20–28° with CF values of 0.28 and 0.25. Manap et al. [21] also mentioned that flat to moderate slopes have higher groundwater probabilities.

In the case of slope aspect, only west (0.22), southwest (0.19) and northwest (0.14) have a positive influence on groundwater occurrence. It is difficult to explain this result, but it might be related to the differences of aspects in terms of vegetation, soil depth and recharge potential. Plan curvature behavior also illustrates that springs are more probable to occur in convex areas. Moghadam et al. [114] also stated that concave areas are more favorable for groundwater compared to the convex slope. Finally, according to the CF value for each influencing factor, karst spring probability map was produced and then categorized into four classes; namely low, medium, high and very high using the natural break classification method (Figure 5), as suggested by Naghibi et al. [31]. Areas with high and very high karst spring potential are seen to be near rivers and faults.

4.2. LR Modelling

According to results of TOL and VIF analysis, no severe multi-collinearity was found between independent variables (Table 2).

As shown in

Table 3. Coefficients of each spring-affecting factors used in logistic regression (LR) model.

Spring-Affecting Factors	β1	S.E2	Wald3	Df4	Sig.5
Geology (Cretaceous)	−12.022	30.011	0.160	1	0.000
Geology (Early Cretaceous)	−3.483	31.228	0.012	1	0.000
Geology (Jurassic-Cretaceous)	8.550	29.001	0.087	1	0.000
Geology (Miocene)	−2.110	33.328	0.004	1	0.000
Geology (Quaternary)	−1.930	30.943	0.004	1	0.000
Geology (Silurian)	10.211	34.550	0.087	1	0.001
Geology (Triassic-Jurassic)	4.685	31.783	0.022	1	0.001
Distance from river (m)	−0.003	0.019	0.025	1	0.000
Slope angle (%)	−0.199	0.021	89.798	1	0.000
Slope aspect (North)	3.340	55.313	0.004	1	0.001
Slope aspect (Northeast)	2.011	51.066	0.002	1	0.015
Slope aspect (East)	4.516	63.014	0.005	1	0.000
Slope aspect (Southeast)	3.782	59.099	0.004	1	0.000
Slope aspect (South)	2.080	58.318	0.001	1	0.001
Slope aspect (Southwest)	13.652	61.021	0.050	1	0.000
Slope aspect (West)	5.318	54.510	0.010	1	0.000
Slope aspect (Northwest)	3.401	50.709	0.004	1	0.001
TWI	−0.990	0.938	1.114	1	0.000
Drainage density (km/km2)	−2.663	0.737	13.056	1	0.000
Altitude (m)	0.374	0.019	387.468	1	0.000
Distance from fault (m)	−0.319	0.044	52.563	1	0.000
Plan curvature	0.044	0.026	2.864	1	0.000
Landuse (Forest)	11.408	73.990	0.024	1	0.012
Landuse (Agriculture)	−0.913	80.428	0.000	1	0.018
Landuse (Rangland)	3.325	72.565	0.002	1	0.011
Constant	−16.33	16.710	0.955	1	0.001

β1= logistic coefficient; S.E.2 = standard error of estimate; Wald3 = Wald chi-square values; Df4 = degree of freedom; Sig.5 = significance.

, results indicated that all applied factors have significant effects on spring occurrence (p < 0.05). Similar results can be found in many previous studies [109].

Table 4. Overall statistics of the LR model with 10 independent variables.

Hosmer and Lemeshow Test		Cox and Snell R2	Nagelkerke R2	Pseudo R2
χ2	Significance
563.31	0.294	0.699	0.821	0.710

showed the overall statistic of the LR model for all 10 independent variables. The results indicated that the goodness-of-fit of the equation could be accepted because the significance of χ² is larger than 0.05 [138]. The values of Cox and Snell R² and Nagelkerke R² showed that the independent variables could explain the dependent variables [138,144].

Using the coefficient values obtained from the final output of the LR analysis, the LR equation was obtained, as shown in Equation (8):

$$\begin{array}{cc}\hfill Z=C_{Geology }+& \left(–0.003\times Distance from river\right)+\left(–0.199\times Slope angle\right)+C_{Slope aspect}\hfill \\ & +\left(–0.990\times TWI\right)+\left(–2.663\times Drainage density\right)\hfill \\ & +\left(0.374\times Topographic elevation\right)+\left(–0.319\times Distance from fault\right)\hfill \\ & +\left(0.044\times Plan curvature\right)+C_{Land use} –16.33\hfill \end{array}$$

(8)

where C_Geology is the LR coefficient for geology factor, whereas C_{Land use} is the LR coefficient for land use factor. All coefficients of Equation (8) were obtained from Table 3. The final generated map using LR model was then classified into low, medium, high and very high classes using the natural break classification method (Figure 6), as suggested by Naghibi et al. [31]. Areas with high and very high potential to karst springs are located in southwest to southeast of the study area.

4.3. Ensemble Modelling

Although LR was capable to implement multivariate statistical analysis, it had some weakness in analyzing the classes of each occurrences conditioning factor. Conversely, the CF model could evaluate the influence of classes of each spring-affecting factor on an occurrence. However, the relationship between karst spring locations and spring-affecting factors was mostly neglected by LR model. As a result, it can be clearly seen that CF and LR models had some weakness when applied individually. The generated map of ensemble CF and LR (ECL) model was prepared and then classified into low, moderate, high and very high classes using the natural break classification method (Figure 7), as suggested by Naghibi et al. [31].

4.4. Validation of Generated KSPMs

The ensemble map was compared to each individual (i.e., CF and LR) map using ROC curve method to find if there was any improvements in combining the individual models. The generated SPMs were overlaid with 30% of spring locations for this phase, which has not been used for model application. According to Figure 8, validation outcomes indicated that LR model could satisfactorily generate groundwater potential map (AUC = 78.4%) and performed better than CF (AUC=66.7%). Moreover, the ECL model had the best performance with an AUC value of 84.7%. It was shown that the ensemble method allowed for a better prediction of groundwater potentiality.

5. Discussion

5.1. Spatial Modeling of Karst Spring Potential

Although several studies have been conducted on construction of karst spring potential maps using statistical models, there are different drawbacks about them, especially defining strict assumptions prior to study. As stated by Dreiseitl and Ohno-Machado [145] and Ozdemir [146], the LR model is usually flexible (especially when it includes only the original covariates) and may be generally desirable. On the other hand, the weak point of bivariate statistical methods such as CF is the neglect of relationship between the groundwater springs and independent variables as observed in this study (i.e., it cannot determine the weight of independent variables). In other word, the CF model just determines the weight of classes of independent variables. Therefore, to overcome these shortcomings, we developed a GIS-based ensemble model that can link the strengths of LR and CF together.

The statistical analysis of LR (Table 3) demonstrated a high contribution of slope aspect (southwest), geology (Cretaceous, Silurian and Jurassic-Cretaceous), land use/land cover (forest) to the modelling process. Prior studies have noted the importance of slope aspect effects on sunshine, air temperature, soil transformation and soil moisture content [147,148]. According to the literature, structural and lithological variations often lead to differences in the permeability and hydraulic properties [8,149,150]. In the case of land use/land cover, it usually influences infiltration and recharge rates as well as the relationship between surface water and groundwater [151]. This is also in agreement with earlier reports, which showed that there is a strong relationship between land use/land cover factor and groundwater recharge [152,153].

5.2. Performance of Ensemble and Individual Models

In Section 4.4, models’ results were analyzed and validated using the validation dataset and the ROC curve method. The results of validation revealed differences between all three approaches, namely including CF (i.e., bivariate statistical model), LR (i.e., multivariate statistical model) and their ensemble.

It should be emphasized that the ensemble model (ECL) provides the best simulation with a good accuracy. This result showed that the accuracy of the prepared KSPM was improved (6.3–17.9%) with the use of the ensemble of bivariate certainty factor and multivariate logistic regression. Therefore, the ensembling process leads to successfully overcome the limitations of LR and CF models. In fact, results shown in this study highlight a remarkable decrease of error through the ensemble approach. This result is in agreement with Arabameri et al. [154] in which the ensemble of LR with evidential belief function (EBF) model can extend an efficient model to reinforce the accuracy of spatial predictions. According to results of this study, the aggregation of multiple prediction of the same variable may lead to better performance, more exact learning and generalization than using a single model prediction. The findings of the current study are consistent with that of Naghibi et al. [30] who found that frequency ratio data mining ensemble model had the best performance among all the single models for groundwater potential mapping. Study in Razandi et al. [26] also compared the performance of different data-mining models and stated that frequency ration (FR) and weight-of-evidence (WOE) could perform better than CF model for groundwater potential mapping in Varamin plain, Iran. This is in line with our finding that the CF model is incapable of analyzing the relation between karst groundwater potential and geo-environmental factors.

6. Conclusions

Springs, as groundwater resources, are important for many sectors, such as domestic consumption and agriculture in arid and semi-arid areas of the world. Sometimes, several families or villages depend heavily on occurrence of a spring; therefore, their spatial modeling is necessary. We designed a new ensemble approach called ECL to deal with this issue. In conclusion, the novel ensemble model gives better predictive performance than could be obtained from any of the constituent models. Therefore, the strength of the proposed methodology is to improve model accuracy by using this novel ensembling technique and to efficiently produce a KSPM for sustainable water recourses management programs. Decision makers and water resources planners can use the proposed approach as an efficient tool to assess the karst spring potential, to establish relationship between karst springs and geo-environmental setting, and to plan and implement effective water resources management in drought condition. Future studies on the current topic are therefore recommended.