Groundwater sustainability assessment based on socio-economic and environmental variables: a simple dynamic indicator-based approach.

The Dehgolan aquifer, which lies in semiarid western Iran, was evaluated using a multi-influencing factor (MIF) analysis to determine groundwater sustainability. Eight indicators, including climatic variability, groundwater exploitation (pumping), groundwater quality, groundwater vulnerability, public participation, legal framework, water productivity, and occupation related to groundwater, were quantified and placed into a series of thematic maps within a GIS framework. Each factor was weighted based on the analyses obtained from the MIF model and the stacked maps were summed to yield a final map showing the degree of sustainability within the groundwater basin. The final groundwater sustainability map showed that 4% of the basin was in a critically unsustainable zone, 30% in an unsustainable zone, 40% in a semisustainable zone, 25% in a sustainable zone, and 1% in an ideally sustainable zone. The final map was validated using a receiver operating characteristic (ROC) method, cross-tabulation, and chi-square tests using groundwater-level decline as a test proxy. The analysis assessed the correlation between water levels that exhibited declines versus the degree of unsustainability of water levels and sustainable water use. The area under the curve was calculated to be 88%, cross-tabulation 64.4%, and the chi-square value was 260.5 with 4 degrees of freedom and values <0.05 (3.627E-55), which suggest that the final map has statistical significance. The sustainability analysis developed is useful as a baseline for development of governance laws to implement management methods in groundwater basins and it can be applied to a wide range of aquifer types in variable climates worldwide.

Introduction

Groundwater is one of the most important natural resources in the world, which today has direct and indirect ties to poverty, food security and economic development, especially in arid and semiarid regions (Konikow and Kendy 2005; Foster et al. 2013; Gleeson et al. 2020; Ostad-Ali-Askari and Shayannejad 2021). Groundwater sustainability has become a significant issue in recent research on world water resources, because more than 2.5 billion people depend on groundwater for drinking water worldwide (UNESCO 2012). Increasing pressure on water resources has led to the rise of new technical research on hydrogeology and sustainable development. ‘Overexploitation’, ‘sustainability’ and ‘vulnerability’ are the most frequently used terms used in literature that assesses the currently stressed groundwater resources (Custodio 2002; Ostad-Ali-Askari et al. 2018, 2019). Perhaps the most common cause of aquifer overexploitation is the lack of adequate or consistent governance in the management of these vital resources resulting in the loss of water supplies. Poor water management also results in creation of environmental stress such as drying of karst springs (Taheri et al. 2015a, 2015c, 2016), salinization of aquifers (Bagheri et al. 2019; Panagopoulos 2021a), brine contamination (Panagopoulos 2020, 2021b), land subsidence (Karimi and Taheri 2010; Rezaei et al. 2021), sinkhole formation (Taheri et al. 2015b, 2019, 2021), or aquifer contamination with chemicals, radioactive materials, and combinations of contaminates that cannot be economically remediated, forcing aquifer abandonment. Unsustainable groundwater management impacts current and future populations and in many cases causes irreparable human hardship.

The desired means of managing groundwater fall within the realm of sustainable development, which has been defined in a variety of ways as discussed in the Brundtland Report (WCED 1987; Brown et al. 1987; Liverman et al. 1988; Dernbach 1998; Parris and Kates 2003; Stoddart 2011; Moore et al. 2017). The simplest and perhaps most quoted of these definitions is the use of resources in a manner that can still meet the needs of the next generation (Harding 2006). A more descriptive definition of sustainable development is “the groundwater extraction regime, measured over a specified time frame, that allows acceptable levels of stress and protects the higher value uses that have a dependency on the water” (Frans et al. 2005). The later definition brings into account the economics of existing use, but also includes a timeframe. Both definitions present concepts that would prevent the destruction of aquifers, such as by permanent compaction causing land subsidence, a reduction in hydraulic conductivity, loss of storage, or permanently changing water quality (e.g., saltwater replacement of freshwater). While overuse may not result in permanent damage to an aquifer, large declines in water levels or changes in water quality can impact local economies for decades. Man-induced impacts on groundwater systems during the Anthropocene have been demonstrated to have impacts similar to past geological events that changed environments, resulting in human migration, mass extinctions and removal of both macro- and micro-ecological systems (Crutzen 2002).

Compaction of porous media and land subsidence are clear signs of unsustainable aquifer management that are easily observed (Khanlari et al. 2012; Xiao et al. 2020; Jiang 2020). However, in many cases, groundwater overexploitation is not so observable. Recognition of unsustainable practices needs to occur before permanent damage to an aquifer system can occur. A first step in evaluating potential harmful practices is to create aquifer vulnerability maps to assess the current state of stress that can be used to develop sustainable management practices. Assessment of current practices can often be contentious within the realm of politics, causing a lack of cooperation between water users and planners. To establish a cooperative environment, socio-economic specialists need to convince water users that unsustainable practices will lead to failure of crops and industries, and could displace populations. Accordingly, designing or presenting evaluation models or multiple integrated maps can provide a comprehensive view of all the factors/agents involved in avoiding a potential crisis or future disaster is very practical and useful. These assessments and maps can enable water managers and planners to educate water users concerning the high potential for disaster occurrence and the magnitude or extent of a current or future crisis.

Numerous research methods have been applied to assess the sustainability of groundwater systems worldwide (Vrba et al. 2006; Gleeson et al. 2010; Pandey et al. 2011; Mays 2013; Hosseini et al. 2019; Majidipour et al. 2021; Samani 2021; Zarei et al. 2022). In many of these studies, the adequacy of the data required to produce a clear view must be evaluated in order to establish if it is sufficient to make the needed assessment (Majidipour et al. 2021). Identifying critical groundwater areas using a variety of factors can provide an overview of sustainability status or a critical overview of an aquifer (Taheri et al. 2020). In addition, the use of classical methods for the assessment of aquifer vulnerability to pollution such as DRASTIC, is a potential first step in making a qualitative groundwater sustainability assessment (Aller et al. 1985; Taheri et al. 2017). Among all of the many methods that have been proposed for groundwater sustainability assessment, the index-based assessment is one of the most objective models that can be used to evaluate the various factors in groundwater sustainability, and can provide a general and understandable view of groundwater status in an area.

Numerous studies have been performed using index-based methods in groundwater management (Bui et al. 2019; Jesiya and Gopinath 2020; Majidipour et al. 2021). Hosseini et al. (2019) reviewed various factors in groundwater sustainability and used 13 of them to evaluate 30 aquifers in Iran. Their database was obtained from the Iran Water Resources Management Company affiliated with the Ministry of Energy (Hosseini et al. 2019). Majidipour et al. (2021) selected 8 different economic and social criteria to study the sustainability of the Mahidasht aquifer and concluded that this aquifer is in a critical state. Perhaps the most universal and practical approach to aquifer assessment is described in the monograph published by Vrba et al. (2006) under the title “Groundwater Resources Sustainability Indicators”. According to the UNESCO report, among the various indicators that could be used for aquifer assessment, 10 specific factors/indicators were suggested for use. These ten factors include renewable groundwater resources per capita, total groundwater extraction/groundwater recharge, total groundwater abstraction/exploitable groundwater resources, groundwater as a percentage of total use of drinking water at the national level, groundwater depletion, total exploitable nonrenewable groundwater resources/annual abstraction of nonrenewable groundwater resources, groundwater vulnerability, groundwater quality, groundwater usability with respect to treatment requirements, and dependence of the agricultural population on groundwater (Vrba et al. 2007). Evaluation of these ten factors can provide an assessment of groundwater sustainability.

A key issue in the evaluation of groundwater sustainability involves the availability of data to support various indicators that can be used in the overall analysis. While several different indicators or factors can be used in a sustainability analysis, solely an increase in the number of factors used does not equate to producing a most robust result in the determination of sustainability (Flour et al. 2021). It is the purpose of this research to use the multi-influencing factor approach to investigate groundwater sustainability in a semiarid region where data are not abundant.

In this study, three groups of sustainability indicators were defined with two to four subfactors placed within each group. A total of eight total indicators were proposed for an assessment of the Dehgolan aquifer in western Iran. A simplistic evaluation was used by applying the multi-influencing factor (MIF). This method has been used by various researchers mainly to evaluate the potential of future groundwater development and/or overexploitation of resources (Magesh et al. 2012; Thapa et al. 2017; Etikala 2019; Nasir et al. 2018; Anbarasu et al. 2020; Etikala et al. 2019; Bhattacharya et al. 2020; Singh et al. 2021; Mandal et al. 2021). Taheri et al. (2020) used this method for the first time in evaluating groundwater critical zones and obtained acceptable results. The first studies using MIF in groundwater investigations were conducted by Shaban (2003) and Shaban et al. (2006). The advantage of this approach is the feasibility of the applied method to obtain a reasonable and acceptable result in arid and semiarid regions with similar conditions to the Kurdistan province study area.

The most important difference between this study compared to the others is the use of a simple and comprehensive methodology so that it is easy to identify the relationship between the eight used indices. In addition, the criteria used, such as employment and climate variability, for the first time in this study, are combined with other factors to assess groundwater sustainability. Nearly all of the equations have been developed by the authors. The other main difference between this study and previous studies (Deines et al. 2021; Samani 2021; Singh and Bhakar 2021) is the compatibility of the model with the minimum available data on various environmental, social and political conditions related to groundwater, so that using the same scores and available data can be applicable for local planning related to groundwater management.

Materials and methods

Study area

The Dehgolan aquifer is the largest, in regards to area, in Kurdistan province, western Iran, which is located to the east of Sanandaj, the capital of Kurdistan province, and within the two counties of Dehgolan and Qorveh. Dehgolan city is the center of the study area and is located 45 km from Sanandaj (Fig. 1a–c). The employment of most people in the county is in agriculture and raising of livestock. The Dehgolan aquifer lies within the southernmost part of Dehgolan-Qorveh and Chardoli (DQCH) subcatchment, which has an area of 2,550 km2 and is a part of the Sefidrud basin. The aquifer system is within the southernmost part of DQCH subcatchments, and commonly is treated as three aquifers—including the Dahgolan, the Qorveh, and the Chardoli aquifers—which are laterally connected. Based on long-term water level observations, hydrographs show that water levels in the aquifers are declining (Fig. 2). Meanwhile, in Chardoli plain, many of the surface effects of groundwater overexploitation are observed, such as land subsidence, resulting in the cracking of buildings, and rise of the well-heads of observation boreholes relative to land surface.

Fig. 1

a Study area in Iran; b the DQCH sub-catchment (D Dehgolan aquifer; Q Qorveh aquifer, and C Chardoli aquifer) and Kurdistan province map; c aquifer border and observation wells; d simplified geological map. Map legend: 1 Quaternary sediments; 2 Hamadan phyllites and Karaj formation; 3 conglomerate and sandstone, gypsiferous marl, shale, flysch turbidite, and calcareous mudstone; 4 Qom formation; 5 Sanandaj shale; and 6 Upper Jurassic granite to diorite (Source: Kurdistan Regional Water Authority)

Fig. 2

Rainfall changes and aquifer hydrographs within the aquifer systems of DQCH subcatchment, including the Dehgolan, Qorveh, and Chardoli aquifers. Groundwater level elevation in m above sea level (m asl)

The study area is the most fertile and mechanized agricultural region in Kurdistan province, which has an important role in the mechanized agricultural economy and agriculture-based–industrial development maintained by increased groundwater pumping. The reduction in water levels within the Dehgolan aquifer, which is the main and largest source of groundwater in Kurdistan, has become a serious challenge for agriculture and industrial development due to groundwater overexploitation and a negative water balance. At least 13,000 ha of farmland are irrigated in this region. Any supply issue in the agricultural sector causes the transfer of uncertainty and instability to the socio-economic systems and industrial development of the region.

The occurrence of cover-collapse sinkholes in recent years in the study area also indicates unsustainability in the groundwater system caused by extraction. If the bedrock was solely limestone, more sinkholes would form. Reports of groundwater contamination in the region with some heavy metals have also been presented in the form of academic theses, but since the results have not been published, they are not cited. All these cases show the importance of this aquifer in the economic development of the region, which extends its effects to other social and environmental aspects.

The Dehgolan plain has an area of 982 km2 including 15% of the DQCH subcatchment. The Dehgolan aquifer has an area of 779 km2 and is the primary source of agricultural water supply for the region. A series of 2,530 deep and semideep (shallower) wells has been installed in the aquifer system. This aquifer is classified as overexploited and is on the list of aquifers that has a ban on new well construction issued by the Ministry of Energy of Iran. This area has an average elevation of 1,876 m above sea level and has a dry and cool climate based on the De Martonne climate classification (Amini and Homayounfar 2017).

Geologically, the study area is located in the Sanandaj-Sirjan structural zone. Location in this zone indicates an abundance of volcanic rocks, metamorphism and the presence of Quaternary magmatism. In the simplified geological map of Fig. 1d, the lithological units of the region are divided into six groups of rocks. These categories are: (1) flood plain, fans and valley terrace deposits, braided and stream channel sediments (Quaternary); (2) tuff, andesitic and basaltic volcanic rocks, tuffaceous shale, slate and meta-limestone and sandstone (Hamadan phyllites and Karaj formation); (3) conglomerate and sandstone; gypsiferous marl, shale, flysch turbidite, and calcareous mudstone; (4) crystalline limestone and dolomite as well as marl, gypsiferous marl, sandy marl and sandstone (Qom formation); (5) dark grey argillaceous shale and schist and phyllite (Sanandaj shale); (6) Upper Jurassic granite and diorite plutonics. However, the oldest rocks in the region are composed of Triassic-Jurassic units (basal units). Quaternary sediments and floodplains are the youngest units in the study area. While some limestone occurs in the area, karst occurrence in the DQCH is not significant. In addition, no detailed studies have been conducted on the limestone and hard rock units in the region and their possible effects on Dehgolan aquifer system groundwater behavior. In the Dehgolan plain, 55 piezometers (observation wells) have been constructed to measure the fluctuations of the groundwater levels, which are measured on a monthly basis (well locations in Fig. 1c).

Indicators selection

Several factors are involved in assessment of groundwater sustainability (GS). Variability in rainfall with a periodicity of 3–5 years can cause temporary lowering of groundwater levels that can complicate a sustainability analysis. Indeed, true climate change can be another factor, which changes the average rainfall in a more permanent manner and impacts a sustainability analysis. Lower rainfall accumulation increases the need for irrigation to maintain crops and puts increasing pressure on aquifer water levels and water quality. However, rainfall alone cannot be the sole factor controlling unsustainability. Anthropogenic water-use practices create additional factors in controlling changing groundwater sustainability. Matching cultivation with the available water resources is critical to groundwater management. Sustainable management of groundwater resources requires adoption of governance to assure that rules and regulations are established and enforced to deal with short- and long-term drought crises and general climate change impacts on aquifer water budgets. Even if adequate rainfall occurs in a basin, other socio-economic factors can encourage excessive water consumption in contradiction to governance laws that oblige sustainable practices.

The entanglement of climate, aquifer hydraulic properties, and socio-economic and governance factors has led researchers to simultaneously study a set of factors involved in GS. Assessment of GS using a set of factors requires a sophisticated approach to help develop governance regulations that are critical in maintaining resources. One approach to this problem is the use of multi-criteria decision models and different data support software within a geographic information system (GIS) framework to evaluate sustainability. The gathering and compiling of the input data for the multi-criteria model is another major issue that must be considered and resolved to produce a proper sustainability analysis. Understanding the aquifer system and water use involves a large number of separate factors to be numerically evaluated. The factors that require evaluation include the hydrogeology of an aquifer system, which includes quantification of the hydraulic properties of unsaturated soils, the saturated zone, the recharge rate, actual exploitation rate (pumping of wells), the permitted discharge of water wells (compliance to rules), identification of unauthorized wells and exploitation control practices, hydraulic relationships between adjacent environments such as karst–alluvium–hard rocks and other geological environments. These are just the technical problems and obstacles in accurately assessing GS programs. When these obstacles are accompanied by socio-economic and political issues, the complexity of the system is greatly increased. Sometimes hasty political decisions cause irreparable damage to the sustainability of aquifers. In Iran, for example, about a decade ago, a law was passed in the Iranian parliament to legalize illegal wells, under which all unlicensed wells were licensed. Political uncertainty cannot be quantified by any model, and all sustainability evaluation models can provide only a shadow of reality.

With all of these limitations in mind, an overall assessment of GS is very useful for local planners and decision-makers because it provides a comprehensive overview. In fact, an assessment of GS is useful to define which factors are well understood, and which factors will require additional study to improve the model, resulting in the development of more realistic governance schemes. As suggested by Vrba et al. (2006) and Hosseini et al. (2019) in past studies, eight indicators have been delegated to evaluate the sustainability of the Dehgolan aquifer. These indicators include three groups and eight indices, which are environmental, socio-political, and economic factors. These properties include climate variability (CV), groundwater exploitation/water-level decline (GEWLD), groundwater quality (GQ), and groundwater vulnerability (GV), socio-political indicators (public participation (PP) and legal frameworks (LF), and economic factors based on occupations related to groundwater use (ORG) and water productivity (WP).

The environmental factors impacted by climate change were adapted from the indicators proposed by UNESCO by Vrba et al.  (2006). Similar socio-political indicators and water productivity were used by Majidipour et al. (2021) to assess GS in a similar environment located in western Iran. The CV and ORG properties are used for the first time in this study as two dynamic factors to assess GS of the Dehgolan aquifer. A conceptual model of the 8 selected indicators in the study is shown in Fig. 3.

Fig. 3

A conceptual model using the eight selected indicators and flow chart in the study

Use of the multi-influencing factor (MIF) method in groundwater sustainability analysis

Multi-influencing factor (MIF) is a mathematical approach used to evaluate the indicators/factors involved in GS assessment. It allows an understanding to be gained on the interactions of the indicators used to assess sustainability. The various factors influencing GS assessment were identified and the interrelationships between these factors were assessed to determine major and minor influences of one factor on another. The MIF method utilizes a weighted overlay analysis, where the weights of different influencing parameters on GS are computed using standard MIF procedures. In the MIF method, every relationship is categorized into major and minor effects, and weighted depending on the impacts of the factor on the other factors. If a major effect exists, showing a strong relationship among parameters, then a score of 1.0 is allocated. If a minor effect exists, revealing weak relationships among parameters, then a score of 0.5 is allocated (Shaban et al. 2006). The relative score for every parameter was calculated from the aggregated weight of both major and minor effects (Table 1). In the next step, the obtained relative score is then used for calculating the proposed score for each selected indicator as shown in Eq. (1).where W indicates the proposed weight of each influencing indicator, A defines the major interrelationship between two indicators, and B is the minor interrelationship between two indicators. (A + B) is the proposed relative rate of effect. Consequently, the relative rates for each influencing factor in GS expressed in scores are as follows:

Table 1

Relative rates and scores for factors influencing groundwater sustainability. Sum values are given in italic

Factor	Major effect (A)	Minor effect (B)	Relative rate(A + B)	% of factor effect or rate:A+B∑A+B×100
Climate variability (CV)	1 + 1	0.5	2.5	13
Groundwater extraction /water-level decline (GEWLD)	1 + 1 + 1 + 1	–	4	21
Groundwater quality (GQ)	–	0.5 + 0.5	1	5
Groundwater vulnerability (GV)	1	0.5	1.5	8
Legal framework (LF)	1 + 1	0.5 + 0.5	3	16
Public participation (PP)	1 + 1	0.5 + 0.5 + 0.5	3.5	18
Water productivity (WP)	1	0.5	1.5	8
Occupations related to groundwater (ORG)	1 + 1	–	2	11
Sum:			19	100

Climate variability (CV): 2 major effects +1 minor = 2(1) + 1(0.5) = 2.5

Groundwater extraction/water-level decline (GEWLD): 4 major effects + no minor = 4(1) + 0 = 4

Groundwater quality (GQ): 2 major effects + no minor = 2(1) + 0 = 2

Groundwater vulnerability (GV): 1 major effect +1 minor = 1(1) + 1(0.5) = 1.5

Legal framework (LF): 2 major effects +2 minor = 2(1) + 2(0.5) = 3

Public participation (PP): 2 major effects +3 minor = 2(1) + 3(0.5) = 3.5

Water productivity (WP): 1 major effect +1 minor = 1(1) + 1(0.5) = 1.5

Occupations related to groundwater (ORG): 2 major effects + no minor = 2(1) + 0 = 2

By using Eq. (1), the percentages of each factor-effect on the GS, after rounding off values, are computed as factor rates. Each indicator is prepared in the form of a spatial map in the GIS environment and each map or layer contains different classes or categories of information (Table 1). Weights of subclasses within each indicator are obtained by the Delphi approach, meaning that they were essentially based on a subjective understanding of the real-world conditions at various well-studied regions. The weights were assumed to be constant, with values as:

w1 = 5 (very high or very influential among other subclasses)

w2 = 4 (high influence)

w3 = 3 (moderate influence)

w4 = 2 (low influence)

w5 = 1 (very low or neutral influence)

Therefore, the factor rates obtained in the MIF method must be multiplied within these category weights. This gives the total weights for the GS final map. All of the indicators with their total weights were incorporated into the GIS environment by applying a weighted overlay analysis in ArcGIS to create the GS maps. These indicators influence weights and subclass rates which were established by the expert panel using a consensus approach. This means that in consultation with the experts of the Kurdistan Regional Water Authority, the priority and latency of these factors and their effect on GS were determined. Weights were then given between 1 and 5 for each subclass as shown in the column heading ‘Score’ in Table 2. The conceptual model of the eight selected indicators is illustrated in Fig. 3.

Table 2

Relative rates and score for factors influencing groundwater sustainability (see text for further detail on the equations)

Factor	Equation	Relative rate	Domain of effect	Qualitative superiority	Score	Total weight
Climate variability (CV) index	CV_Index=Pi-Pav∑i=1n=7Pi-Pav×100	13	–26 to –14	VL	1	13
			–14 to –5	L	2	26
			–5 to 5	M	3	39
			5–20	H	4	52
			20–36	VH	5	65
Groundwater extraction /water-level decline (GEWLD)	GEWLD=GEGE_max×100+WLDWLD_max×100   GEWLD_index=GWLDGWLDmax×100	21	0.3–19	VH	5	105
			19–32	H	4	84
			32–44	M	3	63
			44–58	L	2	42
			58–100	VL	1	21
Groundwater quality (GQ) index	GQI= ECi-ECaECi-ECa×100	5	–99.9 to –30	VH	5	25
			–30 to –5	H	4	20
			–5 to 0	M	3	15
			0–30	L	2	10
			30–100	VL	1	5
Groundwater vulnerability (GV) index	I_DRASTIC=∑i=17Wi×Ri	8	65–75	VH	5	40
			75–85	H	4	32
			85–95	M	3	24
			95–105	L	2	16
			105–117	VL	1	8
Legal framework (LF)	f^x,y=1nh^2∑i=1nKd_ih	16	0–0.2	VH	5	80
			0.2–0.78	H	4	64
			0.78–1.37	M	3	48
			1.37–2.56	L	2	32
			2.56–3.75	VL	1	16
Public participation (PP)	f^x,y=1nh^2∑i=1nKd_ih	18	0–0.1	VH	5	90
			0.1–1.76	H	4	72
			1.76–2.76	M	3	54
			2.76–3.73	L	2	36
			3.73–5.91	VL	1	18
Water productivity (WP) index, by extraction (E)	EWP=In/CWU WPI=Ac×100∑AcEWP×100EWPav	8	0–0.04	VL	1	8
			0.04–0.07	L	2	16
			0.07–0.09	M	3	24
			0.09–0.12	H	4	32
			0.12–0.42	VH	5	40
Occupations related to groundwater (ORG) index	ORGI=NumberofpeoplewhoseoccupationisagricultureTotalpopulationofthevillage×100	11	0–13.5	VH	5	55
			13.5–19	H	4	44
			19–24	M	3	33
			24–41	L	2	22
			41–96	VL	1	11

Validation of an MIF analysis

The economic health of groundwater-dependent jobs is related directly to the groundwater sustainability in a primary used aquifer. Wherein an aquifer is being used unsustainably, there are a set of behavioral issues that permeate society and can disrupt social order. Lack of water tends to increase the rate of rural migration to the cities as shallow and deep wells fail. Therefore, economic health and associated social wellbeing are good criteria for evaluating and validating the MIF in a groundwater sustainability model. The practices of constructing numerous replacement wells, deepening them to the bedrock, and/or replacing pumps with ones having greater lift capacity do not change the underlying issue of water use at unsustainable rates.

Aquifer water level declines can be considered as a proxy in groundwater sustainability assessment in areas where no accurate pumping data are available. This is the case in this study and water levels are used to assess the validity of the developed MIF model. Data from observation wells in the two index years, based on 2008 as the start time during which the aquifer conditions were better than today, and in 2018 when water levels exhibited a deteriorating condition, were used in the assessment. Groundwater level decline rates (WLD) were obtained by calculating the differences in water levels in observation wells between 2008 and 2018. The calculated water-level declines were plotted on the final map at the well locations where measurements were made.

The receiver operating characteristic (ROC) curve method was used to assess the validity of the MIF analysis that was based on a variety of variables. The ROC method has been applied to a wide variety of environmental datasets and is described in detail by Taheri et al. (2021). ROC curve analysis provides a method to differentiate two classes, established through a diagnostic test. This analysis is an assessment of the final distribution of a classification that differentiates between correct and failed predictions according to a 2 × 2 contingency table (Table 3). The confusion matrix is a very useful tool for understanding the results for testing of a model. This matrix shows the results of actual versus predicted class values (Provost et al. 1998). Based on the confusion matrix, a binary target variable creates a “yes” or “no” decision with four possible answers that can occur when assigning a category to the target variable via classification (Wendler and Gröttrup 2021):

Table 3

A confusion matrix for the groundwater sustainability index

		Predicted category
		Yes	No
Actual category	Yes	True positive (TP)	False negative (FN)
	No	False positive (FP)	True negative (TN)

True positive (TP). The true value is “yes” and the classifier predicts “yes”. At a point the aquifer is sustainable and sustainability is computed.

True negative (TN). The true value is “no” and the classifier predicts “no”. A point in the aquifer is unsustainable and no sustainability is computed.

False positive (FP). The true value is “no” and the classifier predicts “yes”. A point within the aquifer is unsustainable, but sustainability is computed.

False negative (FN). The true value is “yes” and the classifier predicts “no”. A point within the aquifer is sustainable, but no sustainability is computed.

These four possible results are typically analyzed within a confusion matrix, which displays the number of TP, TN, FP, and FN. TP, TN, FP, and FN are the counts of true positives, true negatives, false positives, and false negatives, respectively, when the test is applied to large samples (Tables 3 and 4).

Table 4

Relevant metrics of a confusion matrix

Relevant metric	Calculation
Accuracy	TP+TNTP+FN+TN+FP
Sensitivity (TPR)	TPTP+FN
Specificity (FPR)	TNTN+FP
Misclassification rate	FN+FPTP+FN+TN+FP
Precision	TPTP+FP
Balanced Accuracy (AUCb)	^TPR+TNR_2

The most important measures calculated from the contingency table are the false positive rate (FPR) and the true positive rate (TPR). These measures are directly related to the specificity and sensitivity of the outcome, which are used to build an ROC curve. The obvious and most common measure for evaluating an outcome is the ratio of correctly identified data points, the accuracy:with P = TP + FN, the number of actual positive (“yes” or 1) values, and N = TN + FP, the number of actual negative (“no” or 0) values.

The sensitivity and specificity are other basic measures indicating the quality of the model and its ability to predict the individual outcomes correctly. The TPR, or sensitivity, is defined as the ratio of positive values (or 1 s) that are correctly predicted by the model from all actual 1 s. The TPR describes the probability that a positive value data record is classified as such. Thus, a model with a high TPR, near 1, has a high chance of predicting a positive value outcome correctly. If the TPR is low or near 0, then the model is inaccurate in finding the class 1 data, and further parameter tuning or data preparation needs to be considered in order to improve the TPR of the model. The TNR, also called specificity, is the same metric as the TPR, but for the negative value class (or 0 s). More precisely, it is defined as the ratio of negative values (or 0 s) that are correctly predicted by the model from all actual 0 s, which describes the probability that a data record labeled as 0 is recognized as 0 by the model (Wendler and Gröttrup 2021).

The most important values calculated from the contingency table are the false positive rate and the true positive rate. These values are directly related to the specificity and sensitivity of a classifier, and are those used to build an ROC curve. The curve is based on Eqs. (3) and (4) which are:where TP, TN, FP, and FN are the counts of true positives, true negatives, false positives, and false negatives, respectively, when the test is applied to large samples.

In a conventional ROC curve, the true positive rate (TPR) is plotted against the false positive rate (FPR), which is one minus true negative rate. The area under the ROC curve or the AUC is the most important statistical parameter produced in the analysis. The AUC is the primary parameter used for model evaluation. The AUC values of ROC curves provide an assessment of the accuracy of model predictions. Typical values of the AUC lie between 0.5 and 1. A higher AUC value indicates better predictive ability. AUC values less than 0.7 indicate poor predictive capability, with higher values indicating moderate (0.7–0.8), good (0.8–0.9), and excellent (0.9–1) predictive capability (Swets 2014). However, when the AUC value lies above 0.5 there is a degree of statistical significance, but not necessarily a good correlation.

Correlation of the thematic maps and groundwater sustainability indicators (GSI)

The cross-tabulation matrix, which is also known as a contingency table, is the central tool for categorical map comparison (Pontius Jr and Cheuk 2006; Dass 2010). The contingency table, which is used in evaluation in various types of research, is a way to present the frequency distribution of two or more variables concurrently. The different values of the variables are considered to be categories, and they are cross-classified (Malhotra 2007). The cross-tabulation matrix compares two categorical variables by showing a table with the classes of one variable as the rows and the classes of the other variable as the columns (Pontius Jr and Cheuk 2006). Class-by-class paired comparison between the marginal row totals and the marginal column totals allows researchers to evaluate how the two maps relate in terms of the quantity of each class (Pontius Jr 2002; Al Saud 2010). In the current study, spatial correlation between the groundwater sustainability indicator (GSI) map versus eight thematic maps, and GSI versus WLD maps, were measured as an area cross-tabulation, with classes of one map being the rows, and the classes of the second map being the columns.

The chi-square test is used to examine independence across two categorical variables or to assess how well a sample fits the distribution of a known population or goodness of fit (Franke et al. 2012). The formula for computing the test statistic is as follows (Snedecor and Cochran 1980, Eq. 5): where n is the number of cells in the table. Oi and Ei are the observed and expected frequencies for each of the possible outcomes, respectively.

The obtained test statistic is compared against a critical value from the chi-square distribution with (r − 1) × (c − 1) degrees of freedom (df) or (3–1) (3–1) = 4 in this study. Chi-square tests were used to test if the predictors are significant or not (Dixon 2009). At a significance level of 0.05, any value of likelihood less than 0.05 would be significant.

Results

Environmental indicators

Four indicators were selected in this category. Among all of the environmental factors that were considered, climatic variability (CV), the ratio of total groundwater extracted to groundwater level decline (GEWLD), groundwater quality (GQ), and groundwater vulnerability (GV) were used. Climate variability has been touted as a critical factor in the world, and even more recently its effects have been compared to global pandemics such as Coronavirus (Schmidt 2021). Climate variability is like a two-sided coin. Sometimes variability is associated with high rainfall. In a high rainfall case, the rainfall is useful if it has a high associated rate of aquifer recharge, but if it is in the form of flash floods, it does not have a positive effect on aquifer recharge. However, in arid and semiarid regions such as Kurdistan province and western Iran, climate change is generally associated with decreasing rainfall. To evaluate this indicator, a rainfall map was prepared based on long-term data (25 years from 1996 to 2020). This map was divided into five classes; −26 to −14; −14 to −5; −5 to 5; 5–20; and 20–36. (Fig. 4a). The positive impact can be more effective in groundwater sustainability and, vice versa, negative variability has an intense impact on groundwater sustainability and its score is represented by negative values. The equation used by authors to evaluate climate variability is shown in Eq. (6).

Fig. 4

a Climate variability indicator map; the map was reclassified based on the MIF rate (b), and area covered by each class of 5 orders in the CV indicator spatial map (c). d GEWLD indicator map; the map was reclassified based on MIF rate (e), and area covered by each class of 5 orders in the GEWLD indicator spatial map (f)

where Pi is average annual rainfall at station i, and Pav is the average of annual rainfall at all stations. |Pi  Pav| is the absolute value of (Pi – Pav).

The ratio of total groundwater extracted to water level decline (GEWLD) can be modified to a form that defines the total groundwater discharge index to the total of water recharge (Vrba et al. 2006). This indicator has been changed because it is difficult to accurately estimate the spatial aquifer recharge and there was not enough data to accurately calculate this indicator. It should be noted that groundwater extraction is not always equivalent to unsustainability. In places where an aquifer has a high recharge rate, water extraction can have little effect on groundwater level drawdowns. To solve this problem in using this indicator, the ratio of groundwater extraction to the groundwater level decline has been measured. A groundwater extraction map was prepared using interpolated data from a water wells data bank for the study area. To make this parameter dimensionless, the groundwater extraction and water level decline layers were divided by their maximum values based on Eqs. (7) and (8). where GE is the groundwater extraction, GEmax is the maximum value of extracted groundwater, and WLD is the groundwater level decline.

Based on the fact that the higher the WLD and GE values, the greater the likelihood of the numerical value of Eq. (3) approaching 100; therefore, the GEWLD layer was divided into 5 classes—0.3–19; 19–32; 32–44; 44–58; 58–100. These ranges convert to the scores in MIF methods, which are 5, 4, 3, 2, and 1, respectively (Fig. 4c).

The groundwater quality (GQ) indicator was adapted from the groundwater sustainability indicators proposed by Vrba et al. (2006). GQI was determined by using the difference over a 10-year period instead of the ratio of areas with contamination problems to the total area of the aquifer as calculated using Eq. (9) (Fig. 5a). In this approach, changes in electrical conductivity (EC) content have been used as a factor to assess water quality degradation or refreshing. The higher and more positive the percentage, the closer the aquifer is to sustainability, and conversely, high negative values indicate an unsustainable aquifer condition. Due to the lack of accurate data about contaminants in the study area, EC changes over a period of 10 years have been used. It is assumed that human or geological effects have caused increases in EC (increased salinity) and are responsible for deteriorating water quality (Eq. 9). where ECa is actual EC measurements made on water samples from the aquifer, and ECi is the initial EC values measured at the same observation well or data collection location point.

Fig. 5

a GQ indicator map; the map was reclassified based on the MIF rate (b), and area covered by each class of 5 orders in the GQ indicator spatial map (c). d GV indicator based on DRASTIC model map; the map was reclassified based on the MIF rate (e), and area covered by each class of 5 orders in the GV indicator spatial map (f)

Groundwater vulnerability (GV) is a simple definition of the susceptibility of an aquifer to become contaminated, which is assessed based on hydrogeological conditions and examined by one or more different factors. Certainly, the greater the potential of an area to be contaminated, the lower the sustainability of the groundwater resource, based on changes in land use, human impacts, and even environmental effects. Among all of the different methods proposed for aquifer types (karst and alluvial; Taheri et al. 2015b), the DRASTIC method (Aller 1985) is the best approach to assess the vulnerability of alluvial aquifers to contamination based on its universal acceptance and reliable results. DRASTIC is also used as a groundwater sustainability indicator as proposed by Vrba et al. (2006). To use this index, seven different layers were involved in the DRASTIC evaluation (Aller 1985; Taheri et al. 2017). A final vulnerability index map was obtained by the DRASTIC method using the various layer maps and applying Eq. (10). Figure 5c shows the final DRASTIC map of the study area. Based on the different classes of this map, the areas with very high vulnerability, having the most impact on groundwater unsustainability, were given a score of 1. The low impact or very low vulnerability areas were given a score of 5; therefore, a zone with high vulnerability implies a low sustainability state (Fig. 5b).where I is the vulnerability index, Wi is the weighting of each parameter and Ri is the corresponding rate.

Socio-political indicators thematic maps

In addition to all of the environmental factors evaluated to assess groundwater sustainability, political and social issues typically also play very important roles. No scenario for groundwater exploitation is sustainable unless it is supported by a legal framework (LF) and regulatory structure that is implemented and adopted by stakeholder participation. In this study, the existence of dense groups of unlicensed wells has been identified as a weakness of the law in groundwater management. The legal framework in this study refers to whether the existing laws have been able to prevent illegal exploitation or not. The presence of illegal water wells is the best criterion for assessing the strengths and weaknesses of the laws in the study area. The higher the density of illegal wells, the more laws and regulations have failed to play an effective role in groundwater sustainability. To evaluate this indicator, the density of unauthorized wells was used as an unsustainability factor. It is assumed that the higher the number of illegal wells per unit area, the higher the groundwater unsustainability within the area. The impacts of illegal wells were defined based on the GIS help definition Point Density (Spatial Analyst) using a magnitude per unit area from point features that fall within a neighborhood around each cell (Eq. 11). Using the spatial data of illegal wells in the research, a density map of these wells was prepared and divided into five classes. The highest density that has the highest effect on groundwater unsustainability is characterized by the lowest score (1). The lowest number of illegal wells has the greatest positive effect on groundwater sustainability and the highest score (5) (Fig. 6a,b).where (x, y) is the density estimate at the location (x, y); n is the number of observations, h is the bandwidth or kernel size, K is the kernel function, and d is the distance between the location (x, y) and the location of the ith observation (Fotheringham et al. 2000).

Fig. 6

a LF indicator map; the map was reclassified based on the MIF rate (b), and area covered by each class of 5 orders in the GQ indicator spatial map (c). d PP indicator point density map (d); the map was reclassified based on the MIF rate (e), and area covered by each class of 5 orders in PP the indicator spatial map (f)

Public participation (PP) is a critical aspect of developing, implementing, and enforcing any groundwater management scheme that has a primary goal of maintaining groundwater sustainability. Consumers and stakeholders, such as farmers, industry, water companies, and local cooperatives, are directly affected by groundwater deficiency. It is critical that all stakeholders participate in decisions or modifications of groundwater regulations that significantly impact all parties. It should be noted that those groups that do not experience economic impacts are not part of this category. Consumer awareness can contribute to their participation in aquifer management and restoration programs as well as their control of groundwater use. Mostert (2003) considers PP as the direct participation of nongovernment stakeholders in decision-making. Based on available data in the study area, use of volumetric flow meters in the management of water-well extraction is the most accessible criterion for assessing public participation in the implementation of a law or act that can be useful for groundwater sustainability plans. For this purpose, the water wells that have installed flow meters were considered to be a PP positive indicator. This part of the PP layer is quantified by use of spatial data of the wells containing flow meters by use of a density map of flow meter usage. The map was divided into five classes. The highest density has the highest positive effect on groundwater sustainability as characterized by the highest score (5). The low number of flow meters impacts groundwater sustainability, which produces a corresponding low score (1)(Fig. 6d,e).

Economic indicators thematic maps

Water productivity and occupations related to groundwater are two of the important economic factors in groundwater sustainability. Water productivity is the ratio of the net benefits from crop, forestry, fishery, livestock and mixed agricultural systems to the amount of water used to produce these benefits (Molden et al. 2010). In its broadest sense, water productivity states how much water is consumed versus how much economic value the crops generate. Conventional studies of the economic value of water define two types of water productivity, one is physical productivity and the other is economic productivity. The evaluation of the economic value of an agricultural product is a more comprehensive approach, which can be used to assess the production (dollars/unit of water used; Molden et al. 1998). Economic water productivity is defined as the value derived per unit of water used. It can further be used to relate water use in agriculture to nutrition, jobs, welfare and the environment. This concept refers to the value of the product material per cubic meter of water consumed as calculated using Eq. (12).

EWP is the benefit per unit or economic water productivity, ‘in’ is the gross income that is obtained from the sale of a product in a growth season, and CWU is the amount of water used to produce 1 kg of the crop.

To evaluate water productivity in a region, one must first obtain the area under cultivation and the data on the dominant crops in the study area. Based on water well data and research compiled by Amini et al. (2021), the four main crops include wheat, barley, alfalfa and potatoes, which collectively cover about 93% of the irrigated land area of the Dehgolan plain. In terms of economic productivity, wheat, barley and potato crops have the highest net economic benefits and the amount of cultivated land area is in good agreement with the overall crop economic productivity. The unit water values and net benefits for these crops are wheat (1.283 price/cubic meter or p/m3 and 5,066 net income/m3), barley (2.09 p/m3 and 10,980 net income/m3), alfalfa (0.645 p/m3 and 5,494 net income/m3) and potatoes (0.31 p/m3 and 5,345 net income/m3).

The types of crops cultivated in the irrigated areas were determined based on the latest results of the National Census of Water Resources in the Kurdistan Province. In the study area, the water productivity index (WPI) was calculated using Eq. (13). The economic productivity of the four main crops irrigated by wells was multiplied by the total area under cultivation and the result was divided by the total area under cultivation of the different crops. This process can be summarized as follows:

where Ac is the area under cultivation of different crops, ΣAc is the sum of the total area under cultivation in the study area, EWP is economic water productivity, and EWPav is the average economic water productivity for all crops.

Based on the results of the WPI analysis, it was divided into five categories by scores ranging from 1 to 5. After the importance factor was applied, the full range in values was 0.00–0.42. Plots of five different ranges are used in Fig. 7a to show the spatial changes. The divisions used were 0–0.04, 0.04–0.07, 0.07–0.09, 0.09–0.12, and 0.12–0.42 (Fig. 7b).

Fig. 7

a WP indicator map; the map was reclassified based on the MIF rate (b), and area covered by each class of 5 orders in the GQ indicator spatial map (c). d ORG indicator point density map; the map was reclassified based on the MIF rate (e), and area covered by each class of 5 orders in the ORG indicator spatial map (f)

Occupations related to groundwater (ORG) is an indicator used to show the economic importance of groundwater in the employment of rural residents in the study area. This index was calculated based on the number of people working in each village in the agricultural sector. Since employment in the agricultural sector means greater irrigated land based on the data of the Statistics Organization of Iran, the number of people working in the agricultural sector of each village was prepared as a spatial map. To analyze this map, the percent of people employed by agriculture in each village was obtained and divided into five different classes using the following equation (Eq. 14):

The highest percentage of the calculated socio-political indicator in the ORG indicator produces the highest impact on groundwater unsustainability, which is characterized by the lowest score of 1. The lowest number of residents related to irrigated agriculture has the greatest effect on groundwater sustainability and it has highest score at 5 (Fig. 7d,e).

Final groundwater sustainability indicator-based map

The final GSI map for the Dehgolan aquifer was calculated for the MIF model using the weighted sum tool for overlaying the eight indicators in GIS as shown in Fig. 8a. Table 2 shows the final weights resulting from minor and major scores and weights of the three groups of GS indicators. For this purpose, after preparing eight thematic layers related to the three groups of indicators, each layer was reclassified in the GIS environment using the Reclassify tool based on ratings from 1 (lowest importance in the GS) to 5 (highest importance in the GS; Table 2). By multiplying the relative rates in the scores of each class, a weighted map was obtained for the eight indicators. After preparing the eight weighted thematic maps, the final layer was prepared based on the MIF method using the weighted sum command in the GIS package (Fig. 8a).

Fig. 8

a The final map of Dehgolan groundwater sustainability indicators (GSI) by the MIF model; b cover area by each class in the final map. c WLD map; d cover area by each class in the WLD map

Ideally, the sum of the relative rates multiplied by the highest score (5) is 500. The sum of the multiples of the lowest order (1) is 100. The final division of the final GSI map was divided according to the two upper and lower limits and in the form of 400–500, 350–400, 300–350, 250–300 and below 250 (Fig. 8b). The zones or order names used for these classes were: ideally sustainable, sustainable, semisustainable, unsustainable, and critically unsustainable zones, respectively. According to the final map, 4% of the study area is located in the critically unsustainable zone, 30% in the unsustainable zone, 40% in the semisustainable zone, 25% sustainable zone, and 1% in the ideally sustainable zone (Fig. 8b).

Validation

The confusion matrix is a very useful tool for understanding the results for testing of a model. This matrix shows the results of actual versus predicted class values which are used to build an ROC curve. The confusion matrix for the Dehgolan aquifer and its relevant metrics is presented in Table 5. To validate the results of the MIF model in the research, a groundwater-level decline (WLD) map was prepared using data from 54 observation wells in the GIS environment. The maximum rate of WLD in this map varies from −84.4 to 0 m. To produce more accurate results in the validation of the final GSI map, 500 random points were used to create a WLD map and the values of groundwater-level declines and numerical values of the GSI final map were extracted based on these points. In the next step, the error points were removed and finally 494 points were used for the MIF final map validation. In this study, a threshold of −20 m of groundwater-level decline was used as a criterion for validation. Determining the threshold for validation depends on the overall condition of the aquifer. Given that most of Iran’s aquifers have significant declines in water levels, the rate of decline and its significance must be determined based on local conditions. Dehgolan aquifer is also one in which the WLD ranges from a few meters to several tens of meters (Fig. 8c).

Table 5

Application of Eq. (6) and annual precipitation data for calculation of the climatic variability indicator in the study area

Station	Annual precipitation (mm/year)	(Pi – Pav)	|Pi– Pav|	Pi-Pav∑i=1n=7Pi-Pav×100
Dehgolan	288	−11	11	−7.77
Bolban Abad	359.1	60.1	60.1	42.45
Salamat Abad	276.45	−22.55	22.55	−15.93
Khosro abad	307.8	8.8	8.8	6.22
Hasan Khan	298.3	−0.7	0.7	−0.49
Naser Abad	262	−37	37	−26.14
Dosar	300.42	1.42	1.42	1
	Average = 299		Sum = 141.57

Due to the fact that the purpose of this study is to evaluate the groundwater sustainability index, based on the final map, numbers above 350 were determined as sustainable and below it as semisustainable and unsustainable points. All points corresponding to these numbers, which were below −20 m of WLD, were considered as the TP points, while the points that were above 350 in the final map, but had a corresponding WLD of more than −20 m (−20 and more negative points), were designated as the FP points. The FN points are points where the WLD was below the threshold and the corresponding points in the final map were less than 350. The TN points are where the WLD was greater than the threshold, and the model correctly predicted them, and in the final GSI map, the corresponding points are less than 350. To obtain the accuracy of the final map (Eq. 2), TP, TN, FP, and FN values were calculated (Table 6). Based on the predicted values and the actual values of water level decline in Dehgolan aquifer, the accuracy, sensitivity and specificity of the final model were calculated as 83, 73, and 78%, respectively. Other relevant metrics are shown in Tables 4, and 5.

Table 6

Results of the confusion matrix for the groundwater sustainability index in the Dehgolan aquifer, western Irana

		Predicted category
		Yes	No	Total
Actual category	Yes	TP: 91	FN: 33	124
	No	FP: 52	TN: 318	370
Total	–	143	351	494

aRelevant metrics for the Dehgolan GSI map: accuracy = 0.83, sensitivity(TPR) = 0.73, specificity(TNR) = 0.78, misclassification rate = 0.17, precision = 0.64, balanced accuracy (AUCb) = 0.76

To build an ROC curve, for the numbers related to the WLD map range –20 to –84.4, the number 0 was assigned, and for other points (7 to > –20), the number 1 was assigned. Finally, by drawing the ROC curve, the AUC was used to test the final MIF model for the Dehgolan aquifer. In this study, an AUC value of 88% was obtained, which is in an acceptable range based on the proposed categories for this method.

Correlation and cross-tabulation

The objective of the study was to determine zones with optimal groundwater sustainability. For this purpose, the categories of ideally sustainable and sustainable are merged together since these both show areas of sustainable aquifer use. Similarly, critically unsustainable and unsustainable categories are merged together as the ‘unsustainable class’ in the final GSI map. The moderate zone without change is depicted as a separate category. Therefore, the final WLD map can be divided into three categories: sustainable (including very low and low groundwater-level decline) and unsustainable (including very high and high groundwater-level decline) and moderate. The other eight maps were divided into three categories by the same method and each was compared with the GSI final map. The results of cross-tabulation for paired comparison between GSI and overall eight thematic maps and the GSI with WLD maps are shown in Tables 5, 6, and 7.

Table 7

Correlation matrix for the results from the final groundwater sustainability index (GSI) map and groundwater-level decline (WLD) map in the Dehgolan aquifer

			WLD			Total
			USUS	SSUS	SUS
GSI	USUS	Count	127	29	4	160
		% of total	25.9%	5.9%	0.8%	32.6%
	SSUS	Count	53	115	41	209
		% of total	10.8%	23.4%	8.4%	42.6%
	SUS	Count	1	47	74	122
		% of total	0.2%	9.6%	15.1%	24.8%
Total		Count	181	191	119	491
		% of total	36.9%	38.9%	24.2%	100.0%

USUS unsustainable zones, SSUS semisustainable zones, SUS sustainable zones

In the chi-square tests, all variables except ORG have chi-square values less than 0.05, indicating that these features are significant. The chi-square test value for eight variables with 4df based on high to poor significance can be arranged as: GEWLD > PP > LF > CV > WP > GQ > GV > ORG. Based on these results, it is argued that these indicators are significantly related to the final GSI map. In this correlation, it can be concluded that the final GSI map is more sensitive to seven indicators with a chi-square scores less than 0.05 and the final GSI map is less sensitive to the ORG by chi-square scores greater than 0.05 with 4df. However, none of these factors can solely act as the final sustainability index, and it is with the overlap on these thematic maps that the final GSI map becomes meaningful. The final score of the GSI map with WLD shows that the model has good validity, as proved by the ROC and cross-tabulation methods.

Discussion

Recently, several models were proposed for index-based groundwater sustainability assessment (Bui et al. 2019; Hosseini et al. 2019; Majidipour et al. 2021; Zarei et al. 2022). In these groundwater studies, index-based evaluations were more accurate and more applicable than one-dimensional (1D) statistical techniques. Data produced by quantitative indices assessments and the resultant maps are more easily translated into policy-relevant information (Hosseini et al. 2019). In the present study, eight indicators were carefully selected to take into account the environmental, socio-political, and economic characteristics of the Dehgolan aquifer of western Iran for the assessment of aquifer sustainability (Table 2). The first indicator, CV, is an environmental factor and the authors were able to evaluate the climate data to produce a digital map (Eq. 6; Table 4). Climate variability is an undeniable indicator in groundwater management. With the exception of deep well-confined aquifers and fossil groundwater, all groundwater resources are affected by climate change; however, this impact is a function of time and can show up quickly or with a long lag time.

The CV map (Fig. 4a) was weighted by score of 13, which is in the fourth-order of final weights of the MIF model (Fig. 4b; Table 2). Due to the lack of a standard model for classifying the CV map, it was classified based on the environmental break classifier in GIS. This classification method has been used in several other regional groundwater studies (Taheri et al. 2015b and 2020). The spatial extent of the CV map within the MIF method shows that the areas covered by the sustainable and unsustainable zones of GSI with respect to climate variability are 14 and 56%, respectively (Fig. 4b). The semisustainable zone, which is an intermediate state of GS, is neither fully sustainable nor unsustainable, but with a slight change of climate conditions it changes to either condition. The zone covers 30% of the study area. As a result, in terms of climate variability, the Dehgolan aquifer is semisustainable and any change in this index can affect other indicators (Fig. 9a). The correlation of this factor with the final GSI map by the cross-tabulation method is 42.4% (Fig. 9b; Table 6). Based on a chi-square test (χ2 = 62.824) with 4df, the CV vs. GSI maps yielded a p value <0.05 (p = 7.3911E-13) which was statistically significant at the 0.05 level (Table 8).

Fig. 9

a Average percent of area covered by each of the five classes in the eight indicators in the final GSI map; b Cross-tabulation results of GSI vs. eight indicators; c ROC curve of the validated GSI map. Photographs: d–e Shanowarah cover-collapse sinkhole in the Dehgolan aquifer, and overexploitation impacts, including upward exposure of a well casing, suggesting land subsidence (f) and land cracks (g). d–e Photos by Damon Jafarzadeh & Kamal Taheri)

Table 8

The chi-square values and p values of eight thematic maps vs. GSI map, and GSI vs. WLD maps

Indicator correlation	Chi-square value	df	p-value
CV vs. GSI	62.824	4	0.000
GEWLD vs. GSI	291.575	4	0.000
GQ vs. GSI	24.864	4	0.000
GV vs. GSI	24.422	4	0.000
LF vs. GSI	78.687	4	0.000
PP vs. GSI	283.530	4	0.000
WP vs. GSI	38.588	4	0.000
ORG vs. GSI	8.158	4	0.086
GSI vs. WLD	260.461	4	0.000

The GEWLD is the second factor in the group of environmental indicators with a weight of 21 (highest score in the MIF model; Table 2). Due to the fact that the volume of water extracted solely cannot express sustainability/unsustainability, this value (GE) was divided based on the rate of WLD (Fig. 8b) and was proposed as a new indicator, so that more water extraction and greater water-level declines can be measured as unsustainable and vice versa (Fig. 4b). According to the GEWLD spatial extent result, 31 and 43% of the region is in the sustainable and unsustainable states, respectively (Fig. 4b). The correlation of this factor with the final GSI map is 21.4% and shows that this indicator solely is not a comprehensive indicator for assessing GSI (Table 7). The chi-square test for GEWLD versus GSI (χ2 = 291.575) with 4df obtained a p value <0.05 (p = 7.1101E-62) which was significant at the 0.05 level (Table 7).

The GQ is an important indicator in GSI. Different types of contaminants could have been considered for this assessment, but because contamination data were not available, EC data were used as described in the previous sections. Changes in water quality due to salinity increases can be both geogenic and anthropogenic. Therefore, the choice of this factor can be two-sided and its increase is a sign of unsustainability. This factor in the MIF method with weight 5 is in the last order of given weights (Table 2). Based on the GQ spatial map, 52 and 20% are considered as sustainable and unsustainable zones, respectively (Fig. 5a,b). The correlation of GQ with the final GSI map is 28% (Table 6). As a result, the GQ is unreliable for the assessment of GSI without integrating it with other relevant indicators. The chi-square test for GQ versus GSI (χ2 = 24.846) with 4df obtained a p value <0.05 (p = 0.000054) which was significant at the 0.05 level (Table 7).

The GV spatial map was obtained by the DRASTIC method (Aller et al. 1985). It is the fourth indicator and the last in the environmental group indices that shows the inherent potential of the aquifer to become polluted. It is a factor that considers the geological/hydraulic nature of the aquifer in GSI. This layer was divided into five categories based on conventional methods in the DRASTIC model (Fig. 5c). The weight of this layer is eight and within the WP index. It is the sixth in rank of the MIF weights (Table 2). Figure 5b shows 68 and 9% of the total study area covered by sustainable and unsustainable zones, respectively. The correlation of this factor with the final GSI map was 24%, so the GV cannot be used solely to assess overall GS assessment in the Dehgolan aquifer (Table 7). The chi-square test for GV versus GSI (χ2 = 24.422) with 4df achieved a p value <0.05 (p = 0.000066), which is significant at the 0.05 level (Table 8).

The LF with a weight of 16 is located in the third order of the MIF model weights (Table 2). It is the first factor in the socio-political indicators. The control of illegal wells in the sustainable management of groundwater is very important, particularly since the extracted water by these wells is not considered in water budget studies. The LF spatial extent showed 77 and 4.5% as sustainable and unsustainable areas, respectively (Fig. 6a,b). The PP indicator is the second index of the socio-political indicators with a weight of 18 in the MIF model (Table 2; Fig. 6c). The correlation of LF and PP indicators with the final GSI map were 23 and 60.7%, respectively (Table 6). The high correlation of the LF and PP indicators with the final GSI map shows that these indicators are very significant in the GSI mapping. The chi-square tests for LF versus GSI (χ2 = 78.687) and PP versus GSI (χ2 = 283.530) with 4df attained p values <0.05 (p = 3.3052E-16, and 3.8632E-60 for LF and PP, respectively), which were significant at the 0.05 level (Table 8).

The WP is the first index in the economic indicators group. This indicator with a score of 8 is an important factor in evaluation of groundwater sustainability, which was evaluated by calculating the economic productivity of the main products of the region based on Eq. (13). There are many limitations in using this index such as the lack of land cadastral maps; however, in this study, with the available data, the spatial layer of this indicator was prepared and divided into five classes. The WP sustainable and unsustainable zones are 13 and 56% of the spatial extent, respectively (Table 7a,b). The ORG is last index from the economic indices, which has a weight of 11 in the MIF model, and it shows the distribution of groundwater-related employees (Table 2). This map was also divided into five classes and was based on the percentage of population with jobs related to groundwater. It varies from the lowest to the highest at 96% (Fig. 7c). As a rule, the higher the percentage, the higher the groundwater unsustainability expected and vice versa. The spatial map of this index shows 52% in sustainable zones and 14% in unsustainable zones (Fig. 7d). The correlation of WP and ORG indicators with the final GSI map were 27.8 and 30% (Table 7). The chi-square test for WP versus GSI (χ2 = 38.588) with 4df attained a p value <0.05 (p = 8.4742E-8), which was significant at the 0.05 level. However, the chi-square value of ORG versus GSI (χ2 = 8.158) with 4df attained a p value >0.05 (p = 0.085974) which was not significant at the 0.05 level for ORG (Table 8), which indicates the complex nature of the ORG factor and it is not possible to assess sustainability simply by relying on groundwater-related employment data. In other words, even if some of the resident population employment is not related to groundwater, water is still consumed. This fact indicates lack of good governance on the groundwater resource in such areas.

Based on final GSI map, only 1% (7.8 km2) of the total area is located within an ideally sustainable zone (score: 400–442), and 25% (194.8 km2) is located in the sustainable zone (score: 350–400). The semisustainable zone covered 40% (311.6 km2) of the Dehgolan aquifer. The unsustainable and critical unsustainable zones occupied 30% (233.7) and 4% (31.1) of the aquifer, respectively (Fig. 8a,b). In using the GSI map, policy-makers need to lump these indicators together to avoid increases in the area of unsustainable water use. The results of the final GSI map in an overall assessment (percentage of sustainable zones) show that the Dehgolan aquifer, with an area of only 26% for sustainable water use, is essentially an unsustainable aquifer (Fig. 9a). The evaluation of each of the eight indicators without overlaying them, and based on the area covered by these indicators, showed that the GV and LF were associated with the sustainable zones, the PP and ORG were semisustainable, and WP, GEWLD, GQ and GV were unsustainable (Figs. 9a and 10).

Fig. 10

Correlation matrix for the results from the final GSI map versus WLD map, and GSI map vs. the eight indicators in the Dehgolan aquifer

Validation of the final GSI map using the WLD rate and use of three statistical methods, including the ROC curve (88%), cross-tabulation (64.4%) and chi-square value 260.461 with 4df and p-values <0.05 (p = 3.627E−55), showed that the final map has high value. Further, it can be used as a valid road map for future aquifer management programs applied to the Dehgolan aquifer, including restoration and water-based development plans. In addition, cross-tabulation and chi-square tests were used to assess the correlation of different factors with the GSI map. The results of cross-tabulation among each indicator versus the GSI map include, from the highest to lowest, PP (60.7%), CV(42.4%), ORG (30%), GQ (28%), WP (27.8), GV (24%), LF (23.6%), and GEWLD (21.4%; Fig. 9b). These results showed that indicators PP and CV have the most significant correlation with the GSI map, and the sensitivity of these factors is greater than the other factors. These results show that none of the factors can be used solely to provide an overview of the sustainability of the aquifer. In chi-square tests, all factors except ORG have a strong correlation with the GSI map and this suggests the indicators were selected properly (Table 8).

Conclusions

It must be acknowledged that there are many uncertainties in the assessment of groundwater sustainability. Some of these uncertainties and ambiguities are caused by the hidden nature of groundwater and the interaction of various factors that sometimes have no apparent relationships but actually affect each other. In this study, knowing all the complexities of Dehgolan alluvial aquifer system in Kurdistan province, Iran, an attempt was made to evaluate the groundwater sustainability status of this region in a simple way. An important issue in GS evaluation is the understanding of the government experts and the opinions of decision-makers, along with acknowledgment of the limitations of the methods used and the complex relationships between model components. Many new methods, despite their usefulness, are virtually incomprehensible to local decision-makers and managers; however, in this study, using eight different indices that are in the three groups of environmental, socio-political and economic indicators, an index-based groundwater sustainability map of Dehgolan aquifer was prepared. The choices within this method were made with the aim of comparing the model results with the realities in the field, in order to suggest an acceptable method for similar areas in a semiarid climate. Indicators were selected based on the various factors/agents that affect local groundwater and its management, so that these data can be prepared based on current conditions, taking into account the data scarcity in the region. A key objective was to collect data that could be quantified and mapped.

The results showed that 34% of the aquifer area has an unsustainable groundwater resource use, 40% is semisustainable, and 26% is sustainable. The semisustainable zone has a temporary element, as it can become an unsustainable zone with the continuation of the current state of illegal wells and overexploitation. The current trend should constitute an alarm for water managers that regulate the Dehgolan aquifer and indicates the unsustainable situation for the aquifer. Validation of output of this model (GSI map), by use of different statistical methods (AUC area in an ROC curve, cross-tabulation, and chi-square tests) with the WLD map and direct observations of cover-collapse sinkholes, land cracks, and subsidence (Fig. 9d) in the unsustainable zone, indicates the high degree of validity of the MIF model and the reliability of its results. These results, along with the dependence of the Dehgolan county residents on agriculture jobs, conclusively show that if the groundwater in the Dehgolan aquifer is not properly managed, social conflicts and unemployment will increase, thereby causing a harmful breakdown in the economy. According to the results of the study, the highest intensity of unsustainability was documented to occur in the center and southeast portions of the aquifer spatial area.

The groundwater balance of the study area is negative, causing declining water levels, and confirming the unsustainable use of the aquifer. The declining water levels of the region as a whole can be used to confirm the unsustainable use of groundwater but cannot provide spatial analysis on the variations in impact. The MIF model does provide a high degree of spatial analysis of the aquifer condition.

In this study, an innovative method using cross-tabulation tests was applied to evaluate the sensitivity of different layers to the final map and their correlation with the final GSI map. The results showed that climate change and public participation impacts have the highest correlation with the GSI map. Other factors are at lower levels. In the chi-square test, all factors except ORG had a strong correlation with the final GSI. The reason for the weak correlation of ORG with the final model is the uncertainty of the nature of agriculture-related employment and the associated social complexities. The occurrence of cover-collapse sinkholes at Shanowrah in the unsustainable zone, and the confirmation with statistical validation, indicates the high validity of the MIF model. The results of this research demonstrate that the use of simple statistical methods, along with the understanding of aquifer conditions by experts, managers and stakeholders, can be effective in managing groundwater resources and in resolving future disputes and ambiguity surrounding groundwater sustainability.