Effective Use of Catastrophe Multicriteria Decision Analysis in Delineating Groundwater Recharge Potential Zones

Multicriteria decision analysis based on Thom’s [1] Catastrophe Theory (CT) has been applied extensively in solving various social, physical and natural sciences problems. It has become a key tool for identifying the groundwater potentiality of an area with a Geographic Information System (GIS). This paper aimed to apply catastrophe theory effectively by utilizing the standardization formulae suggested by Li et al. [2] and the modified formulae provided by Hongwei Zhu [3]. Depending on the nature of a hydrogeological parameter, CT formulae were chosen to ensure that the least important features also get attention during standardization using Hongwei Zhu's formulae where applicable, which were not possible using Li’s formulae. The standardized values of the features of all the themes obtained using four formulae were sorted in ascending order to estimate the normalized values to ensure that no normalized value of a feature exceeds the others having very close but lower standardized values. The effective use of CT technique was employed in a GIS environment to delineate groundwater recharge potential zones of the middle-west part of Kushtia district, Bangladesh, by integrating influential recharge factors, such as land type, slope, drainage density, distance from surface water bodies, soil permeability, surface clay thickness, rainfall, topographic ruggedness index, topographic curvature index, topographic wetness index and topographic position index. The groundwater recharge potential zones were classified as very good (28.76%), good (32.17%), and moderately good (39.05%) for effective CT technique. But in case of the improper use of CT covering area were 17.24%, 54.05% and 28.71% respectively, and the respective most sensitive factors are rainfall and drainage density. Finally, recharge potential zones were validated using groundwater recharge data with a determinant coefficient of 0.92 and 0.84 for effective and improper use of CT technique respectively.

Introduction

In decision-making, different ratings/weights are set based on importance in a system. Improper assigning of ratings may result in unrealistic and unreliable outcomes. Weighting methods are mainly subjective and objective. The decision-makers assign weights for criteria based on their knowledge and personal judgment in subjective methods. Generally, subjective methods, namely Weighted Linear Combination (WLC), Analytical Hierarchy Process (AHP), Weighted Aggregation Method (WAM) and Weighted Sum Model (WSM), are the most widely used in decision-making. In contrast, the objective method, such as Catastrophe Theory (CT), employs mathematical models to calculate weights that are not impacted by any bias introduced by the decision-maker.

Most natural changes are constant, gradient, continuous or smooth, but they go through a complex process and can be explained and solved using calculus. However, some changes are sudden and transitional, like landslides, earthquakes, volcanic eruptions, etc. A jump-like and discontinuous shift class is called ‘phase-transitions’ or ‘catastrophes.’ Calculus fails to describe or explain these sudden phenomena. Catastrophe theory, first proposed by Thom [1], offers a practical approach to examining phenomena such as sudden changes, multiple modes, and hysteresis. The theory deals with the ‘state variable’ and the ‘control variable.’ CT proposes a unique mathematical relation between the state and control variables. This theory explains how slight and gradual modifications in control parameters can lead to abrupt and discontinuous impacts on the state variables or the dependent variable.

Yang et al. showed that CT significantly reduces the level of subjectivity involved in decision-making [2]. The increasing application has made CT a well-established method in decision-making, extending to the field of decision-making for natural resources, particularly water resource development and management, in recent years. It has been widely applied in various water resources studies, such as the detection of hydrological data discontinuities [3], flood risk evaluation [4], quality of water [5], accessing water security [6], groundwater potentiality by different authors [7]–[12].

CT is used in the catastrophe method to calculate Catastrophe Fuzzy Membership Functions. These functions are obtained by normalizing bifurcation sets. Standardization is crucial in the initial stages of the catastrophe method, as input factors often have varying units of measurement. Raw data must be standardized to a dimensionless quantity ranging from 0 to 1 to enable comparison among different thematic layers. Proposed standardization formulas [13] have been widely adopted by many authors, including [14], [15].

The standardized values of the features of a theme lie between 0 and 1. In most cases, the zero standardized value of a feature does not reflect the actual importance of that feature in that theme and overall performance. In the groundwater potential assessment study [16], the feature index value of 102.33 of the rainfall theme had a standardized value of 0. But no matter how small it is, rainfall influences groundwater potential, but it was neglected. Similarly, authors considered the standardized values of 1 and 0 for the index values of 0.9 and 0.705 for the recharge features [17], which were ill-treated. Similar disparities were noticed in many other studies as well. This indicates that these two formulae are not sufficient for calculating a realistic standardized value from the index of a theme. So, the themes whose features should not have zero standardized value may be standardized using the formulae [18] to overcome the unrealistic ranking.

Another problem arises when standardized values are sorted in descending order. The normalized value of the lower-valued features of a theme produces more influence than a higher-valued feature for closer standardized values [19]. Therefore, the organization of standardized values needs to follow a specific rule before getting normalized values of the features of a theme. The normalization formulae [20] have been utilized in Catastrophe theory depending on the number of control and state variables. Many studies used these normalized formulae [21]–[23]. Authors calculated the normalized values after sorting standardized values in descending orders irrespective of the formulae used for standardization [17], [19]. In contrast, authors [21], [22] did not sort standardized values of different features [13]. This problem of greater influence of the lower-valued feature can be minimized if standardized values are sorted in ascending order before using any catastrophe model for normalization. This problem can be solved by assigning the highest power to a control variable with the lowest value and the lowest power to the control variable with the highest value.

This study proposed an effective version of the catastrophe method by (i) utilizing the standardization formulae [13] and [18] considering the influence of applicable features of a theme and (ii) estimating the normalization value of a feature of a theme using CT models by sorting the standardized values in ascending order to demarcate groundwater recharge potentiality of the middle-west part of Kushtia district, Bangladesh considering as a case study area.

Study Area

The research approach was employed in the central-western region of the Kushtia district, situated on the northern side of the southwest part of the country. The investigation area encompassed three Upazilas (Sub-districts), specifically, Bheramara, Mirpur, and adjacent parts of Daulatpur in Kushtia district, extended to 539.82 km², as depicted in Fig. 1. It is located between 23° 45′ 08″–24° 07′ 52″ N and 88° 51′ 53″–89° 06′ 18″ E geographical coordinates and comprises various villages and two Upazila towns.

Fig. 1. Location map of the case study.

The Ganges (Padma) river and its distributary, the Hisna, are the area’s predominant surface water sources. The Ganges River flows through the north-northeast side of the study area, while the Hisna River entered from the west boundary to the central part and then directed towards the south. The region predominantly comprises deltaic silt, with the northern portion consisting of alluvial sand deposits and the southern part of deltaic sands. The study area is not industrialized and relies on groundwater-based agriculture due to reduced surface water flow in the Ganges River and its distributaries. Thus, the area’s development primarily depends on the agricultural system’s groundwater exploration.

Materials and Methods

Various data sources are utilized in this study, including SRTM imagery, borehole lithology records, groundwater monitoring, and soil texture data. The following information presents the types of data and their origins:

1. The Geological Survey of Bangladesh (GSB) provided Geological Survey data;
2. SRTM image (30 m × 30 m resolution) from USGS;
3. The groundwater monitoring, lithology, and pumping test data were obtained from the Bangladesh Agricultural Development Corporation (BADC), Kushtia;
4. Soil permeability and texture data from Bangladesh County Almanac (BCA) database;
5. Landsat-7 imagery captured on February 6, 2010.

These data were processed, analyzed, and interpreted quantitatively using ArcGIS10.4.1 to evaluate the groundwater recharge potential spatially. The thematic maps of land type (LT), slope (SL), drainage density (DD), surface water body (SWB), soil permeability (SP), surface clay thickness (SCT), rainfall (RF), topographic ruggedness index (TRI), topographic curvature index (TCI), topographic wetness index (TWI) and topographic position index (TPI) were integrated to obtain groundwater recharge potential index (GWRPI). Finally, the area was categorized to find different groundwater potential zones based on GWRPI values.

Preparation of Thematic Layers

Thematic layers of LT, SL, DD, TRI, TCI, TWI, and TPI were generated with the SRTM data in ArcGIS environment. For DD, TR, TWI, and TPI calculations, (1)–(4) have been used. where T_d represents the drainage length; A is the area under study. where E_mean, E_min and E_max are the processing and neighborhoods cells’ mean, minimum and maximum elevation. where α refers to the contributing area per unit contour length in an upslope direction, while β represents the slope angle at the local scale. where M₀ is the elevation of the processing cell, M_n is the grid elevation, n is the total neighborhood cells.

(1) $$\begin{array}{r}DD=T_d/A\end{array}$$

(2) $$\begin{array}{r}TRI=(E_{mean}-E_{min})/(E_{max}-E_{min})\end{array}$$

(3) $$\begin{array}{r}TWI=\ln \left(\frac{\alpha }{tan\beta }\right)\end{array}$$

(4) $$\begin{array}{r}TPI=M_0-\sum _0^{n-1}\frac{M_n}{n}\end{array}$$

Distance from Surface Water Body

Distance from surface water bodies can be utilized to identify the presence of groundwater in an area. The normalized difference vegetation index (NDVI) has been developed as a novel approach for creating maps of surface water bodies. The NDVI utilized visible light (red) and reflected near-infrared radiation to highlight surface water bodies and eliminate soil and terrestrial vegetation features:

(5) $$\begin{array}{r}NDVI=(NIR-VIS)/(NIR+VIS)\end{array}$$

where NIR and VIS represent the near-infrared and visible (red) bands, respectively. Surface water bodies reflect more visible red than near-infrared, resulting in negative NDVI values for water bodies [24].

Surface Clay Thickness

The top geologic formation has a pronounced influence on subsurface water. This study assessed the surface clay-silt-sand thickness using 37 lithological logs. These logs confirmed the presence of clay mixed with silt and sand. A geospatial map of this layer was then created by interpolating the discrete lithological point data.

Rainfall

The replenishment of groundwater depends on the quantity and rainfall duration [25]. There is no rainfall station in the area studied. However, the area’s rainfall was assessed by analyzing data from 35 stations across Bangladesh. The average annual rainfall for 23 years (2000–2022) was interpolated to create a spatial map of the area’s rainfall.

Soil Permeability

Soil is a critical parameter in defining potential recharge to aquifers. The coarse-grained soil can easily infiltrate the surface water due to its high permeability, whereas low-permeable fine-grained soils limit infiltration. This study identified the topsoil textures from the Bangladesh BCA database, whereas the permeability values of these textures are taken from book [26].

Groundwater Recharge Potential Index Estimation

GIS utilizes hydrogeological conditions as the primary mapping units for a given area to create maps of groundwater recharge potential zones. Fig. 2 depicts the methodology employed in delineating groundwater recharge potential zones within the area studied using an effective CT technique. Thematic maps affecting groundwater storage were generated first and then converted into polygons with ArcGIS 10.4.1 software. The union tool combined these thematic maps into a single integrated map. For a comprehensive assessment, the current study utilized a weighted linear combination approach to merge all of the eleven thematic layers, as shown by the following equation:

Fig. 2. Flowchart of groundwater recharge potential zones mapping.

where i^th layer’s weight and rating, denoted by W_j and R_j, respectively, are used in the Equation to determine the GWRPI values for classifying groundwater recharge potential in the area where n represents the number of layers in total. The weight and rating of each layer are critical in the integrated analysis as they significantly impact the final result. The effective CT method was applied to assess and assign weights to all the thematic layers produced.

(6) $$\begin{array}{r}GWRPI=\sum _{j=1}^n{W}_j\times R_j/\sum _{j=1}^n{W}_j\end{array}$$

Multi-dimensional catastrophe fuzzy membership functions can be utilized to resolve inconsistencies in the initial data from various sources [16], [27]. These functions assign values between 0 and 1. The state variable is represented by x, whereas the control variables are denoted by a, b, c and d for the seven catastrophe models shown in Table I.

Table I. Catastrophe Models and Normalization Formulas for CT [20]

Catastrophe model	Control parameter	State variable	Potential function	Normalization formula
Fold	1	1	$V_a\left(x\right)=0.333x^3+ax$	$x_a=a^{1/2}$
Cusp	2	1	$V_{ab}\left(x\right)=0.25x^4+0.5a{x}^2+bx$	$x_a=a^{1/2},x_b=b^{1/3}$
Swallowtail	3	1	$V_{abc}\left(x\right)=0.2x^5+0.333a{x}^3+0.5b{x}^2+cx$	$x_a=a^{1/2},x_b=b^{1/3}$ and $x_c=c^{1/4}$
Butterfly	4	1	$V_{abcd}\left(x\right)=0.17x^6+0.25a{x}^4+0.33b{x}^3+0.5c{x}^2+dx$	$x_a=a^{1/2},x_b=b^{1/3},x_c=c^{1/4}and{x}_d=d^{1/5}$
Oval umb. point	3	2	$V_{abc}(x,y)=x^3+y^3+axy+bx+cy$
Elliptic umb. point	3	2	$V_{abc}(x,y)=x^3-x{y}^2+a(x^2+y^2)+bx+cy$
Parabolic umb. point	4	2	$V_{abcd}(x,y)=x^{2}y+x^4+a{x}^2+b{y}^2+cx+dy$

The four commonly used catastrophe models, namely fold, cusp, swallowtail, and butterfly, are presented in Table I. They demonstrate their normalization forms, which account for all potential discontinuities governed by up to four variables. Three main steps must be implemented in the CT.

Indicator Selection

The selection of factors for the causative sub-system in this investigation depends on data availability. The system and sub-system components are closely interconnected. For instance, an area with lower steepness is more suitable for percolation to groundwater than areas with higher steepness.

For evaluating the groundwater recharge potentiality of the study area, groundwater is treated as a system, whereas LT, SL, DD, SWB, SP, SCT, RF, TRI, TCI, TWI and TPI are considered sub-systems.

Standardization

Standardization is pivotal in converting data into dimensionless values, facilitating comparison and ranking. To standardize the data, the following formulae were used [13]:

For smaller, the better:

(7) $$\begin{array}{r}{x}_j^{{ }^{\prime }}=(x_{j(max)}-x_j)/(x_{j(max)}-x_{j(min)})\end{array}$$

For larger, the better: where symbol j represents the index, whereas x_j signifies the initial value associated with the index j. Additionally, x_j(max) and x_j(min) correspond to each feature class’s upper and lower bounds, respectively.

(8) $$\begin{array}{r}{x}_j^{{ }^{\prime }}=(x_j-x_{j(min)})/(x_{j(max)}-x_{j(min)})\end{array}$$

For the case where the lower index is more significant but the highest value can’t be reduced to zero [18]:

(9) $$\begin{array}{r}{x}_j^{{ }^{\prime }}=x_{j(min)}/x_j\end{array}$$

In the same fashion, when the higher index is more significant, but the lowest value can’t be reduced to zero [18]:

(10) $$\begin{array}{r}{x}_j^{{ }^{\prime }}=x_j/x_{j(max)}\end{array}$$

No matter which formulae are used to find standardized values, the feature classes should be rearranged so that standardized values appear in ascending order.

Normalization

Two principles guide the normalization process: complementary and non-complementary [27]. The approach employed in this study involved using the complementary principle to determine how each control variable would evolve with the occurrence of a catastrophe. Here, a, b, c, and d are the sub-system’s control variables in ascending orders, such as a < b < c < d. As per this principle, control variables counterbalance mutually, ultimately converging towards the mean value of x = (x_a + x_b + x_c + x_d)/4 [22]. Conversely, the control variables under the non-complementary rule cannot counterbalance mutually, and the minimum value of the system’s state variable is determined as x = min (x_a, x_b, x_c, x_d) [21]. Table I showcases the catastrophe models utilized for data normalization in this study.

Feature Classifications of Themes

The Jenks classification technique offers the benefit of recognizing genuine classes in the information and partitioning it into a limited number of classes, which minimizes bias [28]. Thus, Jenks optimization is the most suitable approach for classifying data based on the inherent variability in catastrophe models.

Sensitivity Analysis

This tool is useful in evaluating the importance of subjectivity in components and determining the impact of the weight assigned to each parameter. This analysis is accomplished using two methods: map-removal and single-parameter techniques.

Map Removal

The map-removal method involves eliminating one or more input themes simultaneously to evaluate the sensitivity of input parameters for a suitability analysis [29], [30] and is calculated by:

(11) $$\begin{array}{r}S=100\times (GWRPI-GWRP{I}^{{ }^{\prime }})/GWRPI\end{array}$$

This Equation compares two groundwater recharge potential indices, one obtained with perturbations and the other without. GWRPI obtained from all parameters represents the groundwater recharge potentiality index without perturbations, whereas GWRPI’ represents the perturbed index calculated using fewer layers.

The map removal technique was utilized to assess the impact of individual parameters on GWRPI values by eliminating one layer at a time and observing the resulting changes in the outcomes.

Groundwater Recharge Potential Validation

This study validates the GWRPI map by comparing groundwater monitoring data from single-well pumping tests. The formula used to calculate the groundwater recharge (S_gw) for a unit area:

where Δh represents the increase in water level from pre-monsoon to post-monsoon season, and S_y represents the specific yield. The study area has a fluctuating groundwater level within its surface layer, which comprises clay mixed with silt and sand. The specific yield of clay, silt, and sand are 0.03, 0.08, and 0.23, respectively [31]. An average value of specific yields (0.1133) is used for this study. The annual average groundwater recharge has been computed for each of the 13 observatory wells in the study area using (12) for a unit area for the period 2000–2022.

(12) $$\begin{array}{r}{S}_{gw}=\mathrm{\Delta }h\times S_y\end{array}$$

Results and Discussions

Thematic Layers

The process of creating thematic maps and allocating feature ratings for each theme is explained in detail below.

Land Type

The investigated area was categorized into four groups based on their elevations relative to the mean sea level: lowland (−2–10 m), medium-lowland (10 m–16 m), medium land (16 m–24 m) and highland (24 m–27 m) as shown in Fig. 3a. Highlands facilitate higher runoff, resulting in lower infiltration, so lower ratings were assigned for the highland areas (Tables II and III).

Fig. 3. Thematic map of (a) land type, (b) slope, (c) drainage density, (d) surface water body, (e) soil permeability, (f) top clay thickness, (g) rainfall, (h) topographic ruggedness index, (i) topographic curvature index, (j) topographic wetness index and (k) topographic position index.

Table II. Ratings and Weights of Various Thematic Maps and their Features Using Effective CT Technique

Sub-system	Indicators (Code)	Range of index value	Average index value	Standardized value	Normalized value	Equation used	Weight
Land type (m) (B1)	High (C1)	24–27	25.5	0.16	0.40	9	0.68
	Medium (C2)	16–24	20	0.20	0.58
	Medium-low (C3)	10–16	13	0.31	0.74
	Low (C4)	−2–10	4	1	1
Slope (degree) (B2)	Medium (C5)	6.27–21.08	13.68	0.06	0.25	9	0.65
	Medium-low (C6)	3.21–6.27	4.74	0.18	0.57
	Low (C7)	1.73–3.21	2.47	0.35	0.77
	Very low (C8)	0–1.73	0.87	1	1
Drainage density in km/km2 (B3)	High (C9)	8.87–27.18	18.03	0.09	0.3	9	0.68
	Medium-high (C10)	4.90–8.87	6.89	0.23	0.61
	Medium (C11)	2.97–4.90	3.94	0.40	0.8
	Low (C12)	0.20–2.97	1.59	1.00	1.00
Distance from SWB in m (B4)	Less fav. (C13)	100–2650	1375	0.00	0.00	7	0.66
	Favorable (C14)	0–100	50	0.96	0.99
	Highly fav. (C15)	0–0	0	1	1
Soil permeability in m/day (B5)	Very poor (C16)	1E−07–1E−03	0.0005001	0.010	0.1	10	0.44
	Poor (C17)	1E−04–1E−03	0.00055	0.011	0.22
	Fair (C18)	1E−04–1E−01	0.05005	1.00	1
Surface clay thickness in m (B6)	High (C19)	19.59–27.43	23.51	0.28	0.53	9	0.78
	Mod. high (C20)	14.78–19.59	17.19	0.38	0.72
	Medium (C21)	9.96–14.78	12.37	0.53	0.85
	Low (C22)	3.05–9.96	6.51	1.00	1
Rainfall (B7)	Low (C23)	1445–1460	1452.50	0.9808	0.9903	9	0.9962
	Medium (C24)	1460–1468	1464.00	0.9885	0.9962
	Mod. high (C25)	1468–1475	1471.50	0.9936	0.9984
	High (C26)	1475–1487	1481.00	1.0000	1.0000
Topographic ruggedness index (B8)	Very high (C27)	3.11–7.40	5.25	0.00	0.00	7	0.65
	High (C28)	1.10–3.11	2.10	0.40	0.73
	Medium (C29)	0.50–1.10	0.80	0.56	0.86
	Low (C30)	−5.91–0.50	−2.71	1.00	1.00
Topographic curvature index (B9)	Low (C31)	−14.87–-3.00	−8.94	0	0	8	0.63
	Medium (C32)	−3.00–0.11	−1.45	0.31	0.67
	High (C33)	0.11–6.93	3.52	0.51	0.85
	Very high (C34)	6.93–23.69	15.31	1	1
Topographic wetness index (B10)	Poor (C35)	3.23–7.87	5.55	0.00	0.00	8	0.63
	Fair (C36)	7.87–9.82	8.85	0.30	0.67
	Medium (C37)	9.82–12.58	11.20	0.51	0.84
	High (C38)	12.58–20.84	16.71	1.00	1.00
Topographic position index (B11)	Poor (C39)	10.50–71.88	41.19	0.00	0.00	7	0.70
	Fair (C40)	0.18–10.50	5.34	0.67	0.88
	Medium (C41)	−1.99–0.18	−0.91	0.79	0.94
	High (C42)	−22.08–−1.99	−12.04	1.00	1.00

Table III. Ratings and Weights of Various Thematic Maps and their Features Using Improper CT Technique

Sub-system	Indicators (Code)	Range of index value	Average index value	Standardized value	Normalized value	Equation used	Weight
Land type (m) (B1)	Low (C1)	−2–10	4	1.00	1.00	7	0.64
	Medium–low (C2)	10–16	13	0.58	0.83
	Medium (C3)	16–24	20	0.26	0.71
	High (C4)	24–27	25.5	0.00	0.00
Slope (degree) (B2)	Very low (C5)	0–1.73	0.87	1.00	1.00	7	0.72
	Low (C6)	1.73–3.21	2.47	0.87	0.96
	Medium-low (C7)	3.21–6.27	4.74	0.70	0.91
	Medium (C8)	6.27–21.08	13.68	0.00	0.00
Drainage density in km/km2 (B3)	Low (C9)	0.20–2.97	1.59	1	1.00	7	0.71
	Medium (C10)	2.97–4.90	3.94	0.86	0.95
	Medium-high (C11)	4.90–8.87	6.89	0.68	0.91
	High (C12)	8.87–27.18	18.03	0.00	0.00
Distance from SWB in m (B4)	Highly fav. (C13)	0–0	0	1	1	7	0.66
	Favorable (C14)	0–100	50	0.96	0.99
	Less fav. (C15)	100–2650	1375	0	0
Soil permeability in m/day (B5)	Very poor (C16)	1E−07–1E−03	0.0005001	0.00	0.00	8	0.37
	Poor (C17)	1E−04–1E−03	0.00055	0.001	0.10
	Fair (C18)	1E−04–1E−01	0.05005	1	1
Surface clay thickness in m (B6)	Low (C19)	3.05–9.96	6.51	1.00	1.00	7	0.66
	Medium (C20)	9.96–14.78	12.37	0.66	0.87
	Mod. high (C21)	14.78–19.59	17.19	0.37	0.78
	High (C22)	19.59–27.43	23.51	0.00	0.00
Rainfall (B7)	Low (C23)	1445–1460	1452.50	0.00	0.00	7	0.66
	Medium (C24)	1460–1468	1464.00	0.40	0.74
	Mod. high (C25)	1468–1475	1471.50	0.67	0.90
	High (C26)	1475–1487	1481.00	1.00	1.00
Topographic ruggedness index (B8)	Low (C27)	−5.91–0.50	−2.71	1.00	1	7	0.65
	Medium (C28)	0.50–1.10	0.80	0.56	0.86
	High (C29)	1.10–3.11	2.10	0.40	0.73
	Very high (C30)	3.11–7.40	5.25	0.00	0.00
Topographic curvature index (B9)	Low (C31)	−14.87–−3.00	−8.94	0	0	8	0.63
	Medium (C32)	−3.00–0.11	−1.45	0.31	0.67
	High (C33)	0.11–6.93	3.52	0.51	0.85
	Very high (C34)	6.93–23.69	15.31	1	1
Topographic wetness index (B10)	Poor (C35)	3.23–7.87	5.55	0.00	0.00	8	0.63
	Fair (C36)	7.87–9.82	8.85	0.30	0.67
	Medium (C37)	9.82–12.58	11.20	0.51	0.84
	High (C38)	12.58–20.84	16.71	1.00	1.00
Topographic position index (B11)	High (C39)	−22.08–−1.99	−12.04	1.00	1.00	7	0.70
	Medium (C40)	−1.99–0.18	−0.91	0.79	0.94
	Fair (C41)	0.18–10.50	5.34	0.67	0.88
	Poor (C42)	10.50–71.88	41.19	0.00	0.00

Slope

The area exhibited a slope ranging from 0° to 21.08° and was divided into four classifications based on their level of steepness: very low steep (0°–1.73°), low steep (1.73°–3.21°), medium-low steep (3.21°–6.27°) and medium steep (6.27°–21.08°). Tables II, III, and Fig. 3b display the evaluations for each characteristic.

Drainage Density

The DD map created for this study was categorized into four classes based on their relative importance to groundwater: low-density network (0.20–2.97 km/km²), medium-density network (2.97–4.90 km/km²), medium-high density network (4.90–8.87 km/km²) and high-density network (8.87–27.18 km/km²). The feature ratings for each category are presented in Tables II and III, while Fig. 3c shows the corresponding map.

Surface Water Body

The SWB map was classified into three categories depending on their proximity to the nearest water body: namely the water body (0 m–0 m), a distance from 0 m to 100 m is considered a buffer zone, and a distance from 100 m to a maximum of 2650 m for the rest of the area investigated. These categories were labeled as highly favorable, favorable, and less favorable, as illustrated in Fig. 3d. Since the storage areas have a greater influence on groundwater than areas farther away, the SWB was assigned a higher rating, as presented in Tables II and III.

Soil Permeability

The soil of the area was classified into three groups, as depicted in Fig. 3e. Soils with relatively higher permeability, determined by their sand content and formation permeability, were given higher weightage. As a result, the sandy loam unit was assigned the highest weightage, while clay with some silty clay had the lowest potential for groundwater recharge. Silty clay with some loam having the potential for groundwater recharge was assigned an intermediate weightage (Tables II and III).

Surface Clay Thickness

The spatial distribution of the surface clay-silt-sand thickness prepared using lithological logs is presented in Fig. 3f. Tables II and III show the feature ratings of this theme based on their groundwater yield capacity. It is classified into four categories. The low thickness category ranges from 3.05 m to 9.96 m, while the high thickness category varies from 19.59 m to 27.43 m, and the medium thicknesses are 9.96–14.78 m and 14.78–19.59 m.

Rainfall

The prepared geospatial map for rainfall (Fig. 3g) classified the area into four categories ranging from 1445 mm to 1460 mm, 1460 mm to 1468 mm, 1468 mm to 1475 mm, and 1475 mm to 1487 mm extending from the north, south, and west to the east. Here, higher importance was given to the higher values and vice-versa (Tables II and III).

Topographic Ruggedness Index

The ruggedness map prepared using (3) is shown in Fig. 3h. The area’s index value of −5.91–7.40 was classified as −5.91–0.50, 0.50–1.10, 1.10–3.11 and 3.11–7.40. The rating calculation gave the highest weightage for low ruggedness value and vice-versa (Tables II and III).

Topographic Curvature Index

The prepared geospatial map (Fig. 3i) classified the study area into four zones with TCI values of –14.87–−3.00, −3.00–0.11, 0.11–6.93 and 6.93–23.69. Higher importance is allotted for the higher TCI value and vice-versa (Tables II and III).

Topographic Wetness Index

The TWI map prepared from the SRTM image using (4) was classified into four categories: 3.23–7.87, 7.87–9.82, 9.82–12.58, and 12.58–20.84 (Fig. 3j). The higher TWI indicates a higher groundwater recharge potential. The feature ratings of this parameter are shown in Tables II and III.

Topographic Position Index

It measures topographic slope positions and is frequently used for landform classification. TPI values of the area were calculated using (5), and the generated map is depicted in Fig. 3k. The TPI varying from −22.08 to 71.88 was classified into four categories such as −22.08–−1.99, −1.99–0.18, 0.18–10.50 and 10.50–71.88. The high rating is given to the low TPI value and vice-versa (Tables II and III).

Feature Rating Calculation by Catastrophe Theory

Equations (7)–(10) were utilized to standardize the average index values of various features, whereas (7) was utilized to standardize the thematic layers of SWB TRI and TPI. Equation (8) is used for TCI and TWI. For LT, SL, DD, and SCT, the feature ratings were calculated using (9), and the feature ratings of the remaining thematic layers were standardized using (10) in order to signify less important features, though the referred authors neglected less influential of applicable themes. In order to avoid disorderness in normalization values for close standardized valued features of a theme, standardized values were arranged in an ascending order, as shown in Table II.

The swallowtail catastrophe model (Table I) was utilized for thematic layers SWB and SP, while the butterfly model was applied to the rest thematic layers (as shown in Tables II and III). Finaly, groundwater recharge potential zones in the study area demarcated using ratings and weights of themes. In this study, as the control variables compensate for each other, the complementary principle was used.

Table III shows the ratings and weights of each feature and theme estimated using improper use of CT method, which shows a clear difference in ratings and weights calculated using proper use of CT method as in Table II.

Delineation of Groundwater Recharge Potential Zones (GWRPZ)

A combination of thematic layers was used to demarcate the GWRPZ in the study area. The relevance of each theme to groundwater recharge availability was reflected in its respective feature ratings and weights. The raster layers were converted to polygons, and the union tool in ArcGIS 10.4.1 was used to combine the themes for integration, as described in (6), where each theme weight was multiplied by its corresponding feature rating.

The GWRPZ in the area studied were categorized based on the GWRPI values into three distinct categories: moderately good (0.46–0.71), good (0.71–0.76), and very good (0.76–0.96), shown in Fig. 4a. The north, northeast and south sides, comprising 28.76% of the area, were classified as having very good potential, while the middle portion, accounting for 32.17% of the area, had good potential. The remaining areas, mainly in the middle-east parts, were considered moderately good potential, covering 39.05% of the study area.

Fig. 4. Spatial distribution of GWRPI (a) effective use, (b) improper use and validation of GWRPZ, (c) effective use, (d) improper use of the study area.

Using the improper CT method, the GWRPZ in the area were categorized based on the GWRPI values into three distinct categories: moderately good (0.33–0.62), good (0.62–0.72), and very good (0.72–0.94) as well as shown in Fig. 4b. The northern half and south side, comprising 17.24% of the area, were classified as having very good potential, while 54.05% of the area distributed throughout the whole area had good potential. The remaining areas, mainly in the middle-east part, were considered moderately good potential, covering 28.71% area.

Sensitivity Analysis

The study aimed to evaluate the impact of subjective factors on the identification of GWRPI by employing sensitivity analysis. Specifically, this study used thematic layers to assess how they influence recharge potential.

Map-Removal Sensitivity Analysis

The left side of Table IV displays how removing one layer, estimated using (11), affects the GWRPI. The results show that removing RF causes the highest variation in the GWRPI (mean variation index: 17.72%) due to the high average normalized weight derived from CT. In addition, the GWRPI appears to be highly sensitive to the removal of TPI and SCT, as those have indices of 11.53% and 10.20%, respectively, their weights being 0.77 and 0.71 (Table II). Removing DD causes a variation of 10.03%, whereas SP has the lowest sensitivity with an index of 3.94%, consistent with its low weight. There was a direct correlation between the weights and the mean variation indices. However, some dissimilarity is observed in the case of TWI and SWB.

Table IV. Statistics of the One-Map Removal Sensitivity Analysis

Effective use of CT			Improper use of CT
Parameter removed	Variation Index (%)		Parameter removed	Variation Index (%)
	Mean	SD		Mean	SD
RF	17.72	1.66	DD	12.99	4.12
TPI	11.53	1.42	SL	12.86	4.42
SCT	10.20	2.29	TRI	11.48	3.64
DD	10.03	2.26	TCI	9.74	2.08
LT	9.83	1.95	RF	9.58	4.82
TRI	9.39	2.75	LT	8.83	4.23
SL	8.57	2.61	SCT	8.37	5.25
TCI	8.34	1.33	TPI	7.03	5.63
TWI	5.68	4.35	TWI	6.33	5.07
SWB	4.77	5.43	SWB	5.32	6.30
SP	3.94	3.03	SP	3.12	3.38

The right side of Table IV, generated from the improper use of CT model, displays how removing one layer estimated using (11) affects the GWRPI. The results show that removing DD causes the highest variation in the GWRPI (mean variation index: 12.99%) due to the high average normalized weight derived from CT. In addition, the GWRPI appears to be highly sensitive to the removal of SL and TRI, as those have indices of 12.86% and 11.48%. Removing TCI, RF, LT, and SCT causes a variation of 9.74% to 8.37%, whereas SP has the lowest sensitivity with an index of 3.12%, consistent with its low weight. In this case, more dissimilarities were observed when the weights and mean variation indices were compared.

Validation of Recharge Potential Zones

In this study, the accuracy of the GWRPZ developed using CT was assessed by utilizing the estimated annual average groundwater recharge obtained from the existing 13 groundwater monitoring wells for 2000–2022. The GWRPI map estimated using effective CT was validated by comparing the annual average recharge of 23 years. Groundwater recharge measures how effectively an aquifer system can receive precipitation and surface water. Therefore, the regions with low recharge should reflect low groundwater recharge potential. The recharge zones were compared as very good, good and moderately good. The very good groundwater recharge potential zone belongs to observatory wells 10, 11, 12, and 13. The good groundwater recharge potential zone fits well no. 1, 2, 3, 8, and 9. The rest of the monitoring wells are in a moderately good groundwater recharge potential zone.

But in case of the improper use of CT, the observatory wells 1, 3, 7, and 12 belong to a very good potential zone, whereas wells 2, 5, 6, 10, and 11 fit with good potential zones. The rest wells fall in the moderately good groundwater potential zones.

The average groundwater recharge values are plotted against the GWRPI of three groundwater recharge potential zones, as shown in Figs. 4c and 4d. Here, the recharge values show a linear relation with GWRPI estimated utilizing effective CT technique and validate the groundwater recharge potential zones with a determinant coefficient of 0.92, while the same figure shows a lesser determinant coefficient (0.84) for improper use of CT technique.

Conclusion

This research addresses the effective CT to identify groundwater recharge potential zones and suggests a solution to overcome the problem with improper use of CT technique. Unfortunately, the original standardization formulae of CT provide values from 0 up to 1 for the features of themes. A novel approach was introduced to delineate groundwater recharge potential zones to overcome this conspicuous technical misuse. The effective CT is achieved by (i) using two standardized formulae introduced by Li and other two formulae by Hongwei Zhu, (ii) rearranging the features of each theme in ascending order irrespective of their contribution to each theme, and (iii) using Catastrophe models depending on the number of control and state variables. The recharge potential zones estimated using effective and improper implementation of CT techniques were validated utilizing groundwater recharge data that resulted in a determinant coefficient of 0.92 and 0.84, respectively, suggesting that the groundwater recharge potential zones delineated by the effective CT technique were plausibly reliable to a greater extent. Hence, the proper technique is not only scientifically sound but also reliable and versatile in evaluating groundwater potential. The groundwater recharge potential zones prepared in the present study offer useful suggestions for the government and relevant agencies to ensure the sustainable utilization of groundwater in the region. Since this integrated approach enhances reliability, accuracy, and versatility in evaluating groundwater recharge potential, this cost-effective technique may be applied to any region with similar groundwater influencing factors.