Integrated hydrogeochemical and spatial assessment of groundwater quality using multivariate statistics and GIS: a case study from northeastern Iran

Abstract

Effective monitoring and assessment of groundwater quality are essential for sustainable resource management in arid regions. This study investigates the hydrogeochemical variability and environmental risks associated with groundwater in the Neyshabur Plain, northeastern Iran, as a semi-arid basin under climatic pressures. To clarify spatial variability, the study distinguishes between groundwater inside the main aquifer, representing extraction zones, and outside areas that mainly serve as recharge regions. This separation improves the interpretation of quality patterns and management strategies. A total of 1137 groundwater samples were analyzed to assess physicochemical parameters, major ion composition, and spatial trends across two hydrogeological domains: inside and outside the main aquifer system. Multivariate statistical methods (principal component analysis (PCA) and hierarchical cluster analysis (HCA)) were applied separately for each zone to identify dominant geochemical processes and classify water types. A customized groundwater quality index (GWQI) was developed based on normalized variable loadings and used to map quality gradients using ordinary kriging in a GIS framework. Hydrochemical facies analysis, ion ratio plots, and spatial GWQI maps revealed that more than 70% of samples inside the aquifer fall into poor or very poor quality classes, driven by evaporative dissolution and over-extraction. In contrast, recharge zones show better quality, dominated by carbonate weathering and shorter residence time. The integrated approach provides a replicable model for groundwater quality monitoring and decision-support in data-scarce arid environments. The findings can inform regional environmental policies by identifying vulnerable zones, supporting zoning-based regulation of groundwater abstraction, and prioritizing recharge interventions for long-term aquifer sustainability.

Introduction

Groundwater resources in arid and semi-arid regions are increasingly exposed to a complex array of environmental threats driven by both natural and human-induced factors (Abdalla et al. 2019; Shuaibu et al. 2025). Prolonged droughts, reduced recharge, and rising temperatures associated with climate change exacerbate the vulnerability of aquifer systems. Simultaneously, intensive groundwater abstraction (often unregulated) leads to declining water tables, land subsidence, and increased mobilization of naturally occurring salts (Nagaraj and Masilamani 2023; Chen et al. 2024). In many dry land areas, these processes contribute to the progressive salinization and deterioration of groundwater quality, posing risks not only to agricultural productivity but also to drinking water security and ecosystem stability (Khan et al. 2018; Izeze et al. 2023; Ariman et al. 2024).

Despite growing recognition of groundwater degradation in arid and semi-arid environments, many existing studies remain limited to descriptive analyses of water chemistry or general classifications of water quality (Tawfeeq et al. 2024; Chen et al. 2024; Mohammad et al. 2025). These approaches, while informative, often fail to integrate spatial patterns, statistical drivers, and environmental risk factors in a coherent framework. In particular, the development and application of monitoring-based models (such as regionally calibrated indices, multivariate classification schemes, or GIS-supported zoning tools) are still underutilized in groundwater quality research. The lack of standardized and spatially explicit assessment methods hinders the ability of decision-makers to identify vulnerable zones, anticipate pollution trends, and implement targeted interventions (Nagaraj and Masilamani 2023; Ariman et al. 2024).

The Neyshabur Plain in northeastern Iran exemplifies the environmental and management challenges faced by groundwater-dependent basins in arid regions. The area is characterized by limited and irregular precipitation, high evapotranspiration rates, and geologically complex formations, including widespread evaporative such as gypsum and halite (Naderianfar et al. 2011; Hosseinsarbazi and Ismaeili 2014a; Zeraati Neyshaburi et al. 2022). These conditions create an inherently vulnerable hydrogeological setting in which groundwater is not only the main source of water for domestic and agricultural use but also increasingly threatened by both natural and anthropogenic stressors (Yazdani and Mansouriyan 2014; Rahmati et al. 2015; Amiri Kaiz Kohani et al. 2021). Over the past decades, excessive pumping from deep and semi-deep wells has led to significant water table decline and quality degradation, especially in central and southern sectors of the plain (Hosseinsarbazi and Ismaeili 2014b; Abtahi et al. 2015; Taherian et al. 2021). The lack of continuous monitoring infrastructure, coupled with poor spatial knowledge of aquifer behavior, has further impeded efforts to assess environmental risks and implement evidence-based management strategies. These challenges make the Neyshabur Plain a critical case study for the development of integrated and spatially explicit groundwater assessment tools.

To address these challenges, this study introduces an integrated and zone-specific framework for groundwater quality assessment that combines hydrogeochemical characterization with spatial environmental monitoring. A key innovation lies in the development of a customized groundwater quality index (GWQI), derived from normalized principal component analysis (PCA) loadings, which enables objective and data-driven classification of water quality across a diverse hydrogeological landscape. Unlike traditional indices such as the Canadian Water Quality Index (CWQI) or the weighted arithmetic water quality index, which use fixed weights and predefined parameter sets (Dao et al. 2020; Mohseni et al. 2024), the GWQI developed in this study is data-driven and derived from PCA, ensuring better adaptation to local hydrogeochemical conditions. Moreover, unlike conventional WQIs that may not differentiate between confined and recharge zones, zone-specific GWQI provides more spatially detailed insights, which is crucial for arid aquifers with high variability. In this study, “inside the aquifer” refers to groundwater samples collected from the main Quaternary unconfined aquifer zone with high storage capacity and active groundwater flow. In contrast, “outside the aquifer” includes peripheral zones characterized by shallow or discontinuous groundwater occurrences, typically in semi-confining Neogene formations or areas with limited recharge and low permeability. This operational distinction is based on geological boundaries, piezometric gradients, and lithological transitions observed in the study area. The geostatistical interpolation (ordinary kriging) was applied to visualize spatial patterns of groundwater quality, producing GIS-based maps that highlight critical zones for policy intervention. The combined approach not only enhances scientific understanding of groundwater processes but also offers a replicable tool for resource managers to monitor environmental risks, prioritize intervention zones, and guide sustainable abstraction and recharge planning.

The aim of this study is to develop an integrated environmental monitoring framework to assess groundwater quality in an arid, hydrogeologically diverse region. Specifically, the research seeks to: (1) evaluate the spatial variability of key hydrochemical parameters across distinct groundwater zones using multivariate statistical techniques; (2) construct a statistically weighted GWQI tailored to local conditions; (3) map spatial distributions of groundwater quality using GIS-based interpolation; and (4) support evidence-based groundwater management by identifying vulnerable areas and informing zoning, recharge planning, and abstraction control strategies.

Study area

Geographical setting

The Neyshabur Plain is situated in northeastern Iran, within the boundaries of Razavi Khorasan Province, covering an area of approximately 2400 km² (Fig. 1). It lies between 35°40′ to 36°30′ N latitude and 58°20′ to 59°30′ E longitude, and forms part of the Central Iranian Plateau (Fig. 1.b). The plain is bordered by the Binalood Mountains to the north and northeast, which serve as the primary recharge area for local groundwater systems. Topographically, the plain features a gentle slope toward the south and southwest, facilitating the flow of groundwater and surface water from upland areas to the central basin. This morphology has resulted in the development of extensive alluvial fans and lowland depositional environments, where most groundwater abstraction occurs. The hydrographic network is dominated by ephemeral rivers and seasonal runoff channels originating from the mountainous areas, which contribute limited recharge during wet seasons. The region experiences an arid to semi-arid climate under the Köppen-Geiger classification, with annual precipitation ranging from 200 to 300 mm, concentrated between November and April. Mean annual temperatures range from 13 °C to 18 °C, with summer highs often exceeding 35 °C and winter lows occasionally below 0 °C (Rahmati et al. 2015; Amiri Kaiz Kohani et al. 2021). High evapotranspiration rates (often over 2000 mm/year) combined with limited rainfall result in a severe water imbalance and elevated risk of groundwater salinization.

Geology

The geological structure of the Neyshabur Plain reflects the region’s tectonic complexity and stratigraphic diversity, which have significant implications for groundwater quality and environmental risk. The plain lies within the Binalood geological zone, part of the complex tectonic setting of northeastern Iran, and includes formations ranging from Precambrian metamorphics to recent Quaternary alluvium (Fig. 1.c, Fig. 1.d) (Hosseinsarbazi and Ismaeili 2014b; Abtahi et al. 2015; Taherian et al. 2021). The central and southern areas of the plain are primarily composed of thick Quaternary alluvial deposits, which form the main unconfined aquifer system. These deposits consist of unconsolidated sand, silt, clay, and gravel with varying degrees of permeability. Toward the periphery, particularly in the west and southwest, Neogene marl and conglomerate units act as semi-confining layers, reducing infiltration and altering flow dynamics. In the southeastern and southern margins of the basin, widespread deposits of gypsum, anhydrite, and halite create naturally high levels of total dissolved solids (TDS) in groundwater. These evaporative formations are highly soluble and contribute significantly to the geogenic salinization processes that affect large parts of the aquifer.

[See PDF for image]

Fig. 1

(a) Location of the study area on Iran, (b) location of the study area on Central Flat Basin of Iran, (c) geological map of the study area, and (d) schematic cross section of the study area.

In contrast, the northern and northeastern recharge zones are dominated by fractured limestone and dolomite formations, which exhibit higher secondary porosity and facilitate vertical infiltration of precipitation and runoff. Volcanic and igneous rocks, including andesite and dacite intrusions, are present in the northeastern highlands. While they have limited direct interaction with the aquifer, they influence regional fracture patterns and fault zones that govern groundwater movement and compartmentalization. Structural features such as NW–SE and NE–SW trending faults play a pivotal role in the hydrogeological behavior of the basin. These faults can either enhance groundwater flow through increased permeability or create hydrological barriers depending on their fill and fault gouge characteristics. Their presence explains the localized anomalies in water chemistry and hydraulic gradients observed across the plain (Rahmati et al. 2015; Abtahi et al. 2015; Amiri Kaiz Kohani et al. 2021).

Hydrology

The surface water resources of the Neyshabur Plain consist primarily of seasonal rivers, ephemeral streams originating from surrounding highlands, and man-made irrigation canals. Major rivers such as the Kal-shour and its tributaries flow from the northern and northeastern mountains toward the central and southern parts of the plain (Naderianfar et al. 2011; Rahmati et al. 2015). These rivers are mainly active during the wet season, particularly in late winter and early spring, due to snowmelt and precipitation. Flow rates vary significantly, ranging from several cubic meters per second during peak periods to nearly dry conditions in the summer (Taherian et al. 2021; Zeraati Neyshaburi et al. 2022). Surface water contributes to aquifer recharge primarily through direct infiltration in streambeds (especially in upstream areas with coarse sediments) and through irrigation return flows. Infiltration is more effective in unlined channels and permeable sediments. Water is also managed via small dams, diversion structures, and traditional Qanats. However, due to climatic variability and overexploitation, surface water availability has declined in recent years, intensifying the region’s reliance on groundwater (Hosseinsarbazi and Ismaeili 2014a; Yazdani and Mansouriyan 2014). The Neyshabur Plain aquifer is an unconfined aquifer predominantly composed of Quaternary alluvial deposits, including gravel, sand, silt, and clay, with coarser materials concentrated in the northern and eastern recharge zones (Taherian et al. 2021; Zeraati Neyshaburi et al. 2022). The underlying bedrock, consisting of low-permeability marl, shale, and consolidated limestone, acts as a basal confining layer and plays a key role in controlling vertical groundwater movement and aquifer storage (Naderianfar et al. 2011; Abtahi et al. 2015). The aquifer’s thickness varies across the plain, reaching up to 250 m in certain regions, and serves as the primary water source for agricultural, domestic, and industrial uses. The aquifer undergoes considerable seasonal and long-term fluctuations in water table levels, primarily due to excessive groundwater withdrawal and reduced recharge rates (Fig. 2).

[See PDF for image]

Fig. 2

(a) Location of the sampling sites, (b) piezometric map of the study area (c) types and values of different groundwater uses (d) various types of groundwater discharge and their volumes (Mm³) in the study area

Groundwater recharge is mainly derived from precipitation, mountain streamflow, and irrigation return flows, particularly along the northern and eastern margins (Hosseinsarbazi and Ismaeili 2014a; Zeraati Neyshaburi et al. 2022). Discharge occurs through pumping wells, springs, and natural seepage toward the southern and central areas of the plain. Agriculture is the dominant consumer, extracting approximately 1095 million cubic meters annually. This is followed by drinking water (32.51 Mm³), other uses (15.87 Mm³), and limited combined uses such as agriculture-industry and domestic-agriculture. Industrial and livestock sectors account for less than 3 Mm³. Wells are the principal extraction infrastructure, yielding 102 × 10⁶ Mm³, far surpassing the contributions from springs (611 × 10⁴ Mm³) and Qanats (99 × 10³ Mm³).

Groundwater generally flows from the elevated recharge zones in the north and northeast toward the lower southern and southwestern parts of the plain, following the natural topographic gradient. This flow pattern is governed by the geometry of the aquifer and the spatial distribution of hydraulic conductivity within the alluvial sediments (Rahmati et al. 2015; Amiri Kaiz Kohani et al. 2021; Taherian et al. 2021). According to hydrogeological reports and time-series data from piezometric wells, the plain has experienced a sustained groundwater level decline of 80–100 cm annually over the past two decades, with cumulative drawdowns exceeding 20 m in some central and southern areas (Hosseinsarbazi and Ismaeili 2014a; Yazdani and Mansouriyan 2014; Taherian et al. 2021). This overextraction has resulted in reduced aquifer storage, declining well yields, and in some cases, irreversible compaction and land subsidence. The ongoing depletion of groundwater poses a serious threat to long-term water security, particularly for agriculture, which accounts for over 90% of total groundwater use.

Materials and methods

A total of 1137 samples were analyzed, including 842 wells located within the aquifer boundary and 295 wells situated outside the aquifer system. The dataset was obtained from the Iranian Water Resources Management Company (IWRMC), a governmental authority responsible for groundwater monitoring and management. While detailed well depth data were not available for all sites, records indicate that the average depth of the study wells generally ranges from 150 to 170 m. The chemical parameters measured in each sample included the concentrations of calcium (Ca²⁺), magnesium (Mg²⁺), sodium (Na⁺), potassium (K⁺), carbonate (CO₃²⁻), bicarbonate (HCO₃⁻), chloride (Cl⁻), sulfate (SO₄²⁻), as well as pH, electrical conductivity (EC), TDS, and total hardness (TH). Outliers were identified using box-plot analysis and standard deviation (SD) criteria. Suspicious values were cross-checked and removed if necessary. Missing values were very limited and were handled using simple mean or median substitution depending on data distribution. Electrical neutrality was evaluated by calculating ion balance errors, and samples exceeding acceptable error thresholds (±%5) were excluded from further interpretation (Abdalla et al. 2019; Chen et al. 2024). All chemical concentrations were converted to consistent units (meq L^− 1), and descriptive statistical analyses were performed to examine data distribution. Parameters such as minimum, maximum, mean, and SD were calculated for all measured variables. Subsequently, graphical analyses were performed using Gibbs, Piper, and Durov diagrams to classify water types and identify dominant geochemical processes (Yang et al. 2022; Catania and Reading 2023; Groeschke et al. 2024). Additionally, various ion ratios were computed for each group to further interpret the sources and evolution of groundwater chemistry (Rao et al. 2022; Groeschke et al. 2024). Multivariate statistical analyses, including PCA and (hierarchical cluster analysis) HCA, were then applied independently to each group (Kanji et al. 2025; Mohammad et al. 2025).

Development of a GWQI

The GWQI was developed using a multivariate statistical approach, following methodologies adopted by Abdelaziz et al. (2020) and Ariman et al. (2024). Key hydrochemical parameters were first identified through PCA, where only variables with strong factor loadings (≥ 0.7) on the first principal component (PC1) were selected as dominant indicators of groundwater quality. The relevance of these parameters was further validated using HCA, which showed consistent grouping patterns based on these variables.

Step 1: Determination of weights (W_i).

The weight of each selected parameter was computed as the ratio of its factor loading to the total of all loadings:

1

where, W_i is the weight of the ith parameter, and L_i is the factor loading of the parameter on PC1.

Step 2: Calculation of sub-indices (Si).

Each parameter’s sub-index was computed for every water sample using the following expression:

2

where C_i​ is the measured concentration of parameter i, and S_i​ is its standard permissible limit, based on WHO (2017) drinking water guidelines (Table 2).

Step 3: Calculation of GWQI.

The final GWQI score for each sample was obtained as the weighted sum of the sub-indices:

3

Based on the final GWQI values, groundwater quality was classified as follows (Table 1).

Table 1. Classification of drinking water quality based on GWQI

Range	Water quality for drinking purposes
< 25	Excellent
25–50	Good
50–75	Moderate
75–100	Poor
> 100	Very poor

Software and analytical framework

In this study, various software tools were employed to conduct statistical, hydrochemical, and spatial analyses. IBM SPSS Statistics 21 was used for initial data cleaning, organization, and basic statistical calculations such as minimum, maximum, mean, SD, and ionic ratios. In addition, it was utilized to perform multivariate statistical analyses, including PCA and HCA, to explore patterns and correlations among the water quality parameters. To visualize hydrogeochemical processes and classify water types, Aqua Chem 2014.2 software was employed for constructing Piper, Gibbs, and Durov diagrams separately for the samples inside and outside the aquifer. In this study, Ordinary kriging was applied using default semivariogram setting in, Arc GIS 10.8.2 software to generate spatial distribution maps of groundwater quality parameters. The method was used for visualization purposes. Moreover, Arc GIS 10.8.2 software was used to generate distribution maps of HCA results, and the GWQI.

Results and discussion

Descriptive statistics analysis

The descriptive analysis of the study groundwater samples inside and outside the Neyshabur aquifer is summarized in Table 2. Ca²⁺ concentrations are generally higher in the wells inside the aquifer, with an average of 4.83 meq L^− 1 compared to 3.23 meq L^− 1 in the wells outside. Mg²⁺also follows a similar pattern, with higher average levels inside the aquifer (5.29 meq L^− 1) than outside (4.06 meq L^− 1). Na⁺ exhibits the most significant difference, with the wells inside the aquifer showing nearly double the concentration (19.09 meq L^− 1) compared to those outside (9.42 meq L^− 1).

Table 2. Summary of analytical results of the groundwater samples

	Inside the aquifer				Outside the aquifer
Statistic Parameter	Min.	Max.	Mean	S.D	Min.	Max.	Mean	S.D	Permissible limit (WHO, 2017)
Ca2+ (meq L− 1)	0.30	40.60	4.83	5.51	0.50	28.00	3.23	2.96	5
Mg2+ (meq L− 1)	0.20	38.70	5.29	4.78	0.10	22.80	4.06	3.19	4.1
Na+ (meq L− 1)	0.10	97.30	19.09	18.95	0.10	82.00	9.42	14.50	8.7
K+ (meq L− 1)	0.00	14.00	0.59	1.57	0.00	9.30	0.19	1.00	0.31
CO32− (meq L− 1)	0.00	1.20	0.13	0.25	0.00	1.40	0.20	0.32	-
HCO3− (meq L− 1)	0.70	11.00	2.89	1.24	1.60	8.80	3.91	1.38	-
Cl− (meq L− 1)	0.30	113.70	18.52	22.13	0.20	91.20	7.21	14.36	7.05
SO42− (meq L− 1)	0.15	47.90	3.64	7.40	0.15	41.30	5.70	6.95	5.21
pH	6.60	8.50	8.12	3.26	7.00	7.50	8.28	3.91	-
EC (µS cm− 1)	65.00	14180.00	2939.17	2790.48	4.40	12230.00	1649.99	1954.95	1000
TDS (mg L− 1)	175.77	9128.00	1863.34	1785.98	2.77	7704.90	1039.28	1231.55	500
TH (mg L− 1 CaCO3)	30.00	3200.00	512.85	482.71	60.00	3160.00	373.62	319.70	300

K⁺ levels are low overall in both groups but are higher inside the aquifer (0.59 meq L^− 1) compared to outside (0.19 meq L^− 1). CO₃²⁻levels are slightly higher outside the aquifer (0.20 meq L^− 1) than inside (0.13 meq L^− 1). HCO₃⁻is the dominant anion outside the aquifer, with an average of 3.91 meq L^− 1, compared to 2.89 meq L^− 1 inside. Cl⁻ concentrations are significantly higher inside the aquifer (18.52 meq L^− 1 vs. 7.21 meq L^− 1). SO₄²⁻, on the other hand, is more prevalent outside the aquifer (5.70 meq L^− 1) than inside (3.64 meq L^− 1). pH values are relatively stable across both groups, with a slightly higher average outside the aquifer (8.28) compared to inside (8.12). EC is much higher inside the aquifer (2939.17 µS cm^− 1) than outside (1649.99 µS cm^− 1). TDS also follow this trend, with mean values of 1863.34 mg L^− 1 inside the aquifer and 1039.28 mg L^− 1 outside. TH is likewise greater inside the aquifer (512.85 mg L^− 1 CaCO₃) compared to outside (373.62 mg L^− 1 CaCO₃), consistent with higher concentrations of Ca²⁺ and Mg²⁺ ions in confined zones. In the groundwater samples from the wells inside the aquifer, the dominant cations based on mean concentrations follow the order Na⁺>Mg²⁺>Ca²⁺>K⁺, indicating that Na⁺ is the most abundant. Among anions, the sequence is Cl⁻>HCO₃⁻>SO₄²⁻>CO₃²⁻, suggesting that Cl⁻ is the prevailing anion. Conversely, in the wells outside the aquifer, the cation sequence remains the same (Na⁺>Mg²⁺>Ca²⁺>K⁺), although the overall concentrations are generally lower. However, the anionic sequence differs and follows the order HCO₃⁻>SO₄²⁻>Cl⁻>CO₃²⁻, highlighting HCO₃⁻ as the dominant anion.

Spatial analysis

Figure 3 summarizes the spatial distribution of the physicochemical parameters for the groundwater samples collected from the Neyshabur Plain. In the study area, Ca²⁺ concentrations peak in the central and southern segments of the alluvial aquifer (zones colored yellow–red), while they remain low (dark green) both north of the recharge fringe and beyond the southern margins. Inside the aquifer, Ca²⁺ rises to moderate–high levels (50–100 meq L^− 1) as infiltrating waters dissolve calcite and dolomite within the carbonate-rich fan deposits. Outside the aquifer to the north, freshly recharged mountain‐front springs exhibit Ca²⁺ < 25 meq L^− 1 due to minimal residence time; to the south, saline brines show Ca²⁺ enrichment only where evaporative concentration brings minor gypsum and carbonate dissolution inward. Thus, Ca²⁺ distribution reflects the dual control of residence time in carbonate lithologies (inside) versus fresh recharge or extreme evaporation (outside). Mg²⁺ mirrors the spatial trends of Ca²⁺ but at lower absolute values (typically 5–30 meq L^− 1 inside the aquifer). The highest Mg²⁺ zones coincide with the central basin’s carbonate platforms, confirming dolomite dissolution.

[See PDF for image]

Fig. 3

Spatial distribution maps of the different physicochemical parameters in the study groundwater samples

Beyond the aquifer boundary, Mg²⁺ falls below 5 meq L^− 1 at the mountain-front recharge (north) and only modestly increases (< 15 meq L^− 1) in the southern saline margins, where evaporative drawdown preferentially concentrates the more mobile Na⁺ and Cl⁻ over Mg²⁺. The relative preservation of Mg²⁺ inside versus outside thus highlights the limited interaction of exterior waters with dolomitic source rocks. Na⁺ exhibits its lowest values (< 20 meq L^− 1) in the northern recharge zone, moderate values (20–100 meq L^− 1) throughout much of the aquifer interior, and extreme enrichment (> 200 meq L^− 1) in the southern margin outside the aquifer. Within the aquifer, Na⁺ accumulates through progressive cation exchange (Ca²⁺–Na⁺ swapping on clay surfaces) and minor evaporative concentration. Outside, intensive irrigation return and high evapotranspiration in the low‐lying layers concentrate Na⁺, creating a distinct saline end‐member beyond the mapped boundary. K⁺ remains low across the plain (< 5 meq L^− 1), with only isolated hotspots (5–15 meq L^− 1) linked to irrigated fields where evapo-concentration locally elevates K⁺. Inside versus outside the aquifer boundary, no systematic gradient appears and K⁺ is chiefly governed by soil‐water interactions, rather than regional hydrogeologic flow paths. Zones of elevated CO₃²⁻ appear in the central and southern portions of the aquifer (typically 50–120 meq L^− 1), tapering off both toward the northern recharge fringe (< 20 meq L^− 1) and in the saline margins outside the southern boundary (< 10 meq L^− 1). Inside the aquifer, CO₃²⁻ accumulates as groundwater equilibrates with dissolved carbonate minerals (calcite and dolomite) under relatively high pH conditions, reflecting prolonged water–rock interaction. Outside to the north, fast‐moving, CO₂charged recharge waters remain under-saturated concerning carbonate and thus carry little CO₃²⁻; to the south, high salinity and low Ca²⁺ concentrations limit CO₃² solubility and promote its depletion through common‐ion effects and secondary precipitation. HCO₃⁻ peaks at recharge fringes (300–400 meq L^− 1 north of the aquifer) where CO₂‐rich soil gas promotes aggressive carbonate weathering, then falls to 150–250 meq L^− 1 within the aquifer as secondary calcite precipitates along flow paths. Outside the southern boundary, HCO₃⁻ remains low (< 100 meq L^− 1) because high salinity and evaporation suppress carbonate dissolution. Thus, HCO₃⁻ delineates active recharge versus evaporative‐concentration zones. SO₄²⁻ concentrations are moderate (50–150 meq L^− 1) throughout most of the basin interior (courtesy of gypsum dissolution and oxidation of pyrite in upstream formations) but drop below 25 meq L^− 1 in the recharge zone. A pronounced SO₄²⁻ hotspot (> 200 meq L^− 1) appears just south of the aquifer. This pattern underscores the role of both mineralogy (gypsum/anhydrite) and flow path length in controlling sulfate. Cl⁻ is a sensitive tracer of salinization, with < 25 meq L^− 1 north of the recharge area, 25–150 meq L^− 1 inside the aquifer, and > 300 meq L^− 1 in the southern parts outside the aquifer. The moderate Cl⁻ inside reflects a mix of fresh recharge and minor evaporate dissolution; beyond the boundary, the dissolution of the peripheral evaporate lenses drives Cl⁻ to its highest values, marking a clear salinity front along the aquifer margin. Groundwater pH ranges from mildly acidic to mildly alkaline across the plain. The highest pH values (7.8–8.2) concentrate in the central aquifer, where carbonate buffering is strongest, while the northern recharge fringe exhibits pH near neutral (6.8–7.2) due to a fresh meteoric input and CO₂ dissolution. In the southern margins outside the boundary, pH often drops slightly (6.5–7.0) as the high ionic strength and evaporative enrichment of Cl⁻ and SO₄²⁻ suppress carbonate-buffering capacity. Thus, pH patterns track the balance between fresh, CO₂rich recharge (lower pH), carbonate mineral dissolution (higher pH inside), and concentrated saline brines with weakened buffering (lowered pH outside the southern margin). Both TDS and EC are lowest (< 500 meq L^− 1; < 800 µS cm^− 1) in the northern recharge area, moderate inside the aquifer (500 meq L^− 1; 800 µS cm^− 1), and highest (> 3000 meq L^− 1; >5000 µS cm^− 1) in the southern of outside the boundary. The gradation reflects cumulative mineral dissolution and ion concentration along groundwater flow paths inside the aquifer, versus extreme evapo-concentration in marginal, poorly drained deposits. TH is hard to very hard (200–500 mg L^− 1 CaCO₃) throughout most of the aquifer, driven by Ca²⁺–Mg²⁺ carbonate dissolution, whereas waters outside in the recharge zone are soft (< 150 mg L^− 1 CaCO₃) and those in the saline south are hard due to mixed saline salts. Hardness thus mirrors carbonate interaction intensity and delineates fresh recharge versus mineral‐enriched flow paths.

Mechanism controlling groundwater chemistry

The different hydrogeochemical processes such as evaporation, water-rock interaction, and precipitation can be recognized when plotting the water samples on the Gibbs diagram (Gibbs 1970; Yang et al. 2022). In this diagram, the plot of the ratio of either (Na⁺+ K⁺)/(Na⁺+K⁺+ Ca²⁺) or Cl⁻/(Cl⁻ + HCO₃⁻) as a function of TDS gives a characteristic that is separated into three fields showing different processes affecting the water composition. Based on the Gibbs diagrams for the groundwater samples from the Neyshabur Plain, the dominant natural process influencing water chemistry in both inside (blue) and outside (red) aquifer zones is evaporation-crystallization (Fig. 4). Most samples plot within the evaporation dominance field, with higher TDS concentrations observed inside the aquifer, suggesting more intense mineralization due to longer residence time and reduced recharge (Fig. 4a).

[See PDF for image]

Fig. 4

Gibbs diagram for the study groundwater samples; (a) inside and (b) outside the aquifer

The semi-arid climate of the region, characterized by low precipitation and high evaporation rates, promotes the accumulation of dissolved ions. In addition, the geological composition of the plain, including carbonate-rich alluvium, silicate-bearing formations, and evaporate deposits, supports significant rock-water interactions, contributing to the observed ion ratios. Samples outside the aquifer display relatively lower TDS values, implying younger water or weaker interaction with host rocks (Fig. 4b).

Hydrogeochemaicl classification and water types

The Piper diagram is a widely used graphical tool in hydrogeochemistry that helps classify groundwater types and identify geochemical processes by plotting the relative concentrations of major cations (Ca²⁺, Mg²⁺, Na⁺, K⁺) and anions (Cl⁻, SO₄²⁻, HCO₃⁻, CO₃²⁻) in two triangular fields, which are then projected into a central diamond to determine the overall water facies (Piper 1944; Appelo and Postma 2005; Catania and Reading 2023). The Piper diagrams for the groundwater samples inside and outside the aquifer of the Neyshabur Plain reveal distinct geochemical characteristics and evolutionary trends shaped primarily by natural geological processes (Fig. 5). In both groups, the cation triangle shows dominance of Na⁺ and K⁺ over Ca²⁺ and Mg²⁺, while the anion triangle is heavily skewed toward Cl⁻ and SO₄²⁻ rather than HCO₃⁻. This pattern is more pronounced in samples outside the aquifer, where Cl⁻ and SO₄²⁻ concentrations are particularly high, indicating a more advanced stage of chemical evolution (Fig. 5b).

[See PDF for image]

Fig. 5

Piper diagram for the study groundwater samples; (a) inside and (b) outside the aquifer

The central diamond-shaped field shows that most water types plot in the Na⁺–Cl⁻ and Na⁺–SO₄²⁻ facies, which are typical of groundwater affected by prolonged water–rock interaction, dissolution of evaporate minerals (e.g., halite and gypsum), and possibly ion exchange processes. Inside the aquifer, while Na⁺–Cl⁻ face still dominates, there is a relatively broader distribution with a few samples showing influence from Ca²⁺ and HCO₃⁻, suggesting some mixing with recharge zones or interaction with carbonate-rich formations (Fig. 5a). The high presence of Na⁺ and Cl⁻ across both zones is consistent with the regional lithology, which includes salt-rich and SO₄²⁻-bearing deposits. The patterns indicate that natural geochemical mechanisms, such as silicate and evaporate weathering, ion exchange, and prolonged residence time, play a key role in shaping the groundwater chemistry in this arid to semi-arid basin.

The hydrogeochemical distribution of the main water types (higher than 10 numbers in any class) depicted in the two pie charts offers a clear distinction between groundwater inside the main aquifer, and waters sampled outside the aquifer system, primarily from recharge or upstream zones (Fig. 6). These differences reflect the geochemical evolution of groundwater as it moves from recharge areas to discharge zones.

In the main aquifer zone, Na⁺-Cl⁻ water type is dominant (227 samples), followed by Na⁺-Mg²⁺-Ca²⁺-Cl⁻ (113), and Na⁺-Mg²⁺-Cl⁻-SO₄²⁻ (51), indicating a high degree of mineralization (Fig. 6.a). These water types suggest longer residence time, significant cation exchange, and evaporative concentration. The aquifer system mainly overlaps with fine-grained Quaternary deposits, clay flats, and poorly sorted conglomerates where groundwater flow is slower, and water-rock interaction is more extensive. The abundance of Na⁺ and Cl⁻ ions implies that ion exchange has occurred (most likely involving clay minerals) and that some areas may be affected by evaporation. In contrast, groundwater outside the aquifer shows a higher frequency of Mg²⁺-Ca²⁺-HCO₃⁻ (34 samples), Na⁺-Mg²⁺-Cl⁻ (30), and Mg²⁺-Ca²⁺-Na⁺-HCO_₃⁻ (25) types (Fig. 6.b).

[See PDF for image]

Fig. 6

Pie diagram of the different water types for the study groundwater samples; (a) inside and (b) outside the aquifer

These zones correspond to recharge areas with alluvial fans, carbonate rocks, and volcanic formations in the foothill regions. The dominance of HCO₃⁻-bearing facies in these areas indicates active recharge from precipitation or surface runoff, where water has undergone limited geochemical evolution. The presence of Ca²⁺ and Mg²⁺ in significant quantities is consistent with carbonate dissolution from limestone and dolomite units, while mixed facies with Na⁺ and Cl⁻ suggest initial stages of silicate weathering or influence of volcanic rock dissolution. The differences between the two domains (inside vs. outside the aquifer) can thus be attributed to variations in residence time, geological substrate, and hydrogeological dynamics. Waters in the recharge zone are younger, less mineralized, and closely tied to local lithology, while those in the aquifer are older, more evolved, and influenced by regional flow paths and geochemical processes such as ion exchange, evaporation, and dissolution of evaporate and silicate minerals. The hydrogeochemical facies identified in the study area align with findings from various arid and semi-arid regions globally, where water-rock interactions and evaporative concentration are common processes influencing groundwater geochemistry (Ravikumar and Somashekar 2017; Ouarani et al. 2020; Benyoussef et al. 2024). However, differences arise when comparing with regions indicating recharge from freshwater sources and limited evaporation (Umamageswari et al. 2019; Patel et al. 2023). For instance, the higher prevalence of SO₄²⁻ in certain regions suggests a greater influence of gypsum dissolution in the above studies compared to the Neyshabur Plain.​ Also, the presence of confined aquifers in the mentioned studies introduces additional complexity in hydrochemical facies distribution, which may not be as prominent in the Neyshabur Plain.

The Durov diagram provides a comparative visualization of the hydrogeochemical evolution of groundwater samples from within and outside the Neyshabur Plain aquifer. In Fig. 7.a, corresponding to internal aquifer samples, the majority of data points cluster in the Na⁺–Cl⁻ and Na⁺–SO₄²⁻ facies fields, indicating high mineralization and advanced water–rock interaction processes, likely due to longer residence time, evaporation, and ion exchange in confined aquifer conditions. The spread of points toward the lower right quadrant also suggests cation exchange, where Ca²⁺ and Mg²⁺ are replaced by Na⁺. In contrast, Fig. 7.b, representing samples outside the aquifer, shows a significant concentration in Mg²⁺–HCO₃⁻ and Ca²⁺–HCO₃⁻ zones, reflecting fresher water influenced by recharge from carbonate-rich formations. This suggests active infiltration, limited water–rock interaction, and shorter groundwater flow paths. The broader spread in red samples also implies a wider range of hydrochemical processes in the recharge zones, including the mixing and dissolution of carbonates. The Durov diagrams affirm the pie chart findings by highlighting clear hydrogeochemical zonation between confined and recharge areas. These differences can be attributed to geological variability (e.g., marls and evaporates in the center vs. carbonates and alluvial fans on the margins), hydraulic gradients, and aquifer confinement, which control the flow dynamics and chemical evolution of the groundwater in the Neyshabur Plain. Although the Piper and Durov diagrams may initially suggest that the wells located inside and outside the aquifer share similar hydrochemical characteristics, a more detailed geochemical analysis reveals important distinctions that align with the hydrogeological settings of each group.

[See PDF for image]

Fig. 7

Durov diagram of the study groundwater samples; (a) inside and (b) outside the aquifer

Wells outside the aquifer, located mainly in recharge and upstream zones, tend to show fresher water signatures with higher proportions of Ca²⁺ and HCO₃⁻, reflecting limited residence time and carbonate weathering. In contrast, samples from inside the aquifer display more evolved chemical compositions dominated by Na⁺–Cl⁻ and Na⁺–SO₄²⁻ facies, consistent with prolonged water–rock interaction, ion exchange, and evaporative dissolution. These trends, while not always sharply contrasted in the diagrams alone, become evident when integrated with geological and hydrodynamic context, thereby reinforcing the relevance of the diagrams in distinguishing between the two hydrogeological domains.

Ionic ratios

The comparative assessment of ten key ionic ratios in the groundwater samples from inside and outside the Neyshabur aquifer provides critical insight into the dominant hydrogeochemical processes and lithological influences across the study area (Table 3; Fig. 8). Inside the aquifer, 69% of the samples exhibited Na⁺/Cl⁻ >1, indicating plagioclase (albite) weathering or cation exchange processes (Fig. 8.a). In contrast, 79% of samples outside the aquifer fell in this category, suggesting a stronger influence from silicate weathering in the unconfined zone (Fig. 8.b). The higher frequency of Na⁺/Cl⁻ < 1 inside the aquifer (29% vs. 15% outside) also points to more significant halite dissolution, likely from deeper evaporative-bearing layers. A parallel study in the KomaduguYobe Basin (Sahel, Nigeria) reported Na⁺/Cl⁻ ratios ranging from 0.42 to 2.75 (mean 1.11), with 55% of shallow samples > 1 (meteoric origin) and the remainder affected by saline intrusion or anthropogenic inputs (Shuaibu et al. 2025). Similarly, in the Essaouira Basin (semiarid Morocco), a Na⁺/Cl⁻ ratio close to unity was interpreted as an active halite dissolution under strong evapotranspiration (Tahiri et al. 2016). Both zones in the study area were dominated by samples with Ca²⁺/Mg²⁺ ≤ 1 (66% inside, 64% outside), suggesting dolomite weathering as the dominant process (Fig. 8.c and 8.d). However, a slightly higher percentage of the outside aquifer samples (10% vs. 6% inside) showed Ca²⁺/Mg²⁺ ≥ 2, reflecting more active silicate weathering or calcite enrichment in the shallow aquifer layers. This mirrors findings from the Changbai Mountain region (China), where hydrochemical facies of HCO₃⁻–Ca²⁺–Mg²⁺ and Ca²⁺–Na⁺–Mg²⁺ exhibited predominantly Ca²⁺/Mg²⁺ < 1, reflecting pervasive dolomite weathering in lacustrine–volcaniclastic aquifers (Li et al. 2022). The Cl⁻/HCO₃⁻ ratio clearly distinguishes evaporative vs. recharge conditions. A majority of the inside aquifer samples (65%) exhibited Cl⁻/HCO₃⁻ >1, pointing to high evaporation or possible saline intrusion (Fig. 8.e). Conversely, 67% of the samples outside the aquifer showed Cl⁻/HCO₃⁻ < 1, highlighting fresher recharge-dominated conditions in the unconfined system (Fig. 8.f).

Table 3. Ionic ratios indicating of possible geochemical processes for the study groundwater samples

Ionic ratio a	Value	Possible geochemical interpretation	Inside the aquifer b	Outside the aquifer	Considerations
Na+/Cl−	< 1	Reverse ion exchange or seawater and brine	29	15	Evaluates source of salinity (halite vs. weathering or exchange)
	= 1	Halite dissolution	2	6
	> 1	Plagioclase (albite) dissolution	69	79
Ca2+/Mg2+	≤ 1	Dolomite dissolution	66	64	Distinguishes between carbonate rock types
	= 1–2	Calcite or gypsum dissolution	28	26
	≥ 2	Silicates dissolution	6	10
Cl−/HCO3−	< 1	Freshwater recharge	34	67	Salinity indicator and recharge status
	= 1	Transitional water	1	1
	> 1	Evaporation or saline water intrusion	65	32
SO42−/Cl−	< 0.5	Halite dissolution	33	4	Indicates sulfate source
	= 1-0.5	Balanced source	1	1
	> 1	Gypsum and anhydrite dissolution	66	95
Ca2++Mg2+/SO42 + HCO3−	< 1	Cation exchange or silicate input	72	74	Evaluates contribution of rock weathering vs. ion exchange
	= 1	Carbonate/sulfate balance	2	4
	> 1	Excess Ca2+-Mg2+ from gypsum or calcite	26	22
Na+/Na++Ca2+	< 0.5	Carbonate weathering	16	41	Indicator of water-rock interaction and exchange processes
	= 0.5	Mixed processes	1	1
	> 0.5	Cation exchange, silicate weathering	83	58
Mg2+/Mg2++Ca2+	< 0.5	Calcite weathering	34	36	Differentiates dolomite from calcite dissolution
	= 0.5	Balanced carbonate source	3	4
	> 0.5	Dolomite weathering	63	60
HCO3−/Cl−	< 1	Evaporation	66	31	Reveals recharge condition
	= 1	Transitional regime	1	1
	> 1	Recharge area of freshwater	33	68
SO42−/HCO3−	< 1	Carbonate dissolution	31	61	Identifies dominance of sulfate-bearing lithologies
	= 1	Carbonate and gypsum dissolution	1	1
	> 1	Gypsum dissolution	68	38
TDS/EC c	< 0.5	Presence of organics or low ions	1	0	Quality check; higher values may indicate organics or highly mineralized water
	= 0.5–0.6	Normal range	3	1
	> 0.6	Non-electrolyte or saline effects	96	99

^a Ions concentrations in meq L^− 1

^b Values in percent

^c TDS in mg L^− 1; EC in µS cm^− 1

[See PDF for image]

Fig. 8

Ionic ratios of the study groundwater samples (blue and red points are the inside and outside the aquifer samples, respectively)

Comparable contrasts occur in semiarid Morocco, where saline intrusion facies prevail near the coast but inland-perched aquifers exhibit stronger bicarbonate signatures under active recharge (Bourjila et al. 2024). These parallels underscore how climatic gradients and aquifer confinement universally govern hydrogeochemical evolution in arid settings. A striking contrast appears in the SO₄²⁻/Cl⁻ ratio: 66% of the inside aquifer samples and 95% of the outside aquifer samples had a ratio > 1, suggesting widespread gypsum/anhydrite dissolution across the region, with especially strong influence in the shallower, more oxidized zones outside the aquifer (Fig. 8.g and 8.h). Lower ratios (< 0.5), typical of halite, were more common inside (33%) than outside. Both inside and outside the aquifer, 72% and 74% of samples respectively had (Ca²⁺ + Mg²⁺)/(HCO₃⁻ + SO₄²⁻) < 1, which suggests cation exchange or silicate weathering processes dominating over CO₃²⁻ / SO₄²⁻ equilibrium (Fig. 8.i and 8.j). The slightly higher number of samples with ratios > 1 inside the aquifer (26% vs. 22%) may indicate localized gypsum or calcite inputs. Likewise, in the central Rif, (Ca²⁺ + Mg²⁺)/(HCO₃⁻ + SO₄²⁻) < 1 in most samples, indicating silicate input or ion exchange over pure carbonate dissolution (Benyoussef et al. 2024). A significantly larger portion of inside aquifer samples (83%) had Na⁺/(Na⁺ + Ca²⁺) > 0.5, compared to 58% outside (Fig. 8.k and 8.l). This indicates more advanced cation exchange processes or greater silicate weathering inside the aquifer, likely due to longer residence time and more evolved groundwater chemistry. Conversely, 41% of outside samples had values < 0.5, reflecting the dominance of fresh carbonate input. This strong cation exchange signature under long residence times is echoed in the Bokoya Massif (central Rif), where reverse exchange (Na⁺ gain, Ca²⁺ loss) shapes Na⁺–Cl⁻ facies in deeper wells (Bouaissa et al. 2021). Dolomite weathering appears dominant in both zones, with 63% of the inside aquifer and 60% of outside aquifer samples showing Mg²⁺/(Mg²⁺ + Ca²⁺) > 0.5 (Fig. 8.m and 8.n). This confirms the presence of dolomitic formations in both settings, though with a slightly stronger dolomite signature inside the aquifer. Similar proportions were reported in Sahelian aquifers under mixed evaporative and recharge regimes (Sidibe et al. 2019). Reinforcing earlier observations, 66% of the inside aquifer samples had HCO₃⁻/Cl⁻ < 1 (suggesting evaporation), while 68% of the outside aquifer samples had ratios > 1, confirming recharge-dominated, fresher water conditions outside the aquifer (Fig. 8.o and 8.p). Gypsum dissolution dominates both zones but is especially prominent inside the aquifer, where 68% of the samples showed SO₄²⁻/HCO₃⁻ >1, compared to 38% outside (Fig. 8.q and 8.r). This suggests deeper interaction with evaporative strata inside the confined aquifer. Outside, a higher percentage of the samples (61%) had ratios < 1, consistent with carbonate dissolution and shallow recharge. The TDS/EC ratio exceeded 0.6 in nearly all samples (96% inside, 99% outside), indicating highly mineralized water across both zones (Fig. 8.s and 8.t). These elevated values suggest the presence of non-electrolyte solutes or high ionic strength due to prolonged water–rock interaction in an arid climate.

HCA

HCA performed on the groundwater samples from both inside and outside the aquifer of the Neyshabur Plain has revealed three distinct hydrogeochemical clusters in each domain (Fig. 9; Table 4). Inside the aquifer, cluster 1 (75% of samples) is the dominant cluster and represents the baseline hydrochemical composition. The Stiff diagram for this group shows balanced major cations (Ca²⁺, Mg²⁺) and anions (Cl⁻, SO₄²⁻), indicating typical geogenic influence via water-rock interaction. The prevalence of this cluster (three-fourths of all the inside samples) suggests a well-buffered, geochemically stable aquifer with little anthropogenic interference. It likely represents deeper wells tapping into confined zones with relatively long groundwater residence time. Cluster 2 (22% of samples) shows a shift in chemical signature, possibly due to localized variation in recharge or lithology. Slight elevations in Na⁺ and Cl⁻ suggest ion exchange. Its spatial distribution likely corresponds to transition zones between fresh and mineralized water bodies, possibly at the margins of the aquifer or near mixing zones. Cluster 3 (3% of samples) is chemically identical in both inside and outside groups, highlighting a common hydrochemical type that transcends aquifer boundaries. This may reflect either a shared recharge origin or similar environmental conditions influencing certain locations inside and outside the aquifer. Chemically, it displays moderate ionic concentrations and is positioned between the red and blue clusters in terms of composition. Outside the aquifer, cluster 1 (80% of samples) mirrors cluster 1 inside in terms of proportional representation and overall ionic pattern. However, the hydrochemical homogeneity is even more pronounced here, possibly due to shallower groundwater, shorter flow paths, and stronger surface influence. These conditions can minimize geochemical evolution, resulting in chemically uniform water influenced by direct recharge and surface interactions. Cluster 2 (16% of samples) is chemically more enriched in Na⁺, Cl⁻, and TDS, suggesting anthropogenic inputs such as fertilizer leaching, return flows from irrigation or contamination from surface activities.

Table 4. Identified clusters of the study groundwater samples based on HCA

Inside the aquifer			Outside the aquifer
	Samples (number)	Samples (%)		Samples (number)	Samples (%)
Cluster 1	635	75	Cluster 1	238	80
Cluster 2	186	22	Cluster 2	48	16
Cluster 3	21	3	Cluster 3	9	4

[See PDF for image]

Fig. 9

Dendrogram and Stiff diagrams of the identified clusters of the study groundwater samples; (a) inside the aquifer; (b) outside the aquifer (the mean concentrations of each cluster were used to draw the Stiff diagram; the cluster shown in green is the same in both groups, therefore, a single chart was plotted for it)

This cluster likely reflects external, human-driven processes in unconfined or semi-confined conditions with limited natural filtration. Cluster 3 (4% of samples) is shared between both domains, acting as a hydrochemical bridge. Its presence outside the aquifer in a similarly limited proportion indicates that certain groundwater bodies outside the confined zone undergo processes or have sources similar to those in the interior (perhaps due to overlapping recharge areas or lateral flow).

The HCA results from the Neyshabur Plain, with three hydrochemical clusters both inside and outside the aquifer, are consistent with studies from India, Egypt, Spain, and Morocco (Ramesh and Elango 2012; El sheikh et al. 2021; De la Hera et al. 2022; El Marnissi et al. 2024). A dominant cluster with geogenic signatures, representing 75–80% of samples, is commonly reported in semi-arid regions such as the Nile Delta and Rajasthan (Gad 2015; Singh et al. 2020). Minor clusters showing elevated Na⁺ and Cl⁻, as seen in Neyshabur’s pink and blue groups, are often associated with anthropogenic impacts like agriculture or wastewater (Khan et al. 2019; El Marnissi et al. 2024). The presence of a shared hydrochemical cluster (green) across both zones in the Neyshabur is a relatively rare finding in international studies. Similar cluster patterns based on salinity and ion composition have been observed in Pakistan and Morocco (Ahmed et al. 2018; El Marnissi et al. 2024). Overall, Neyshabur’s results align with global trends but also highlight localized features such as aquifer connectivity and mixed recharge sources.

Figure 10 summarizes the spatial distribution of the clusters for the study groundwater samples. Inside the aquifer, cluster 1 dominates the central and western parts. These areas are generally characterized by thicker alluvial deposits with well-developed groundwater flow systems. This cluster likely reflects mature groundwater with longer residence times, allowing extensive water-rock interactions. Cluster 2 is more scattered and occurs in intermediate zones where geological transitions, such as changes in sediment grain size or mineralogy, may influence groundwater chemistry.

[See PDF for image]

Fig. 10

Spatial distribution of the identified clusters of the study groundwater samples

Cluster 3 appears in localized areas with geochemical characteristics similar to those found outside the aquifer, suggesting zones of hydrogeological continuity or shared recharge conditions. Outside the aquifer, cluster 1 is mostly located in upstream and peripheral regions. These zones are often dominated by younger, less evolved groundwater derived from direct infiltration of runoff and precipitation, with minimal residence time. Cluster 2 represents areas where local lithology, such as the presence of marls, conglomerates, or fine-grained sediments, affects the composition of the groundwater. Cluster 3, shared with the aquifer interior, indicates that in some areas, natural recharge and geological conditions lead to similar groundwater evolution on both sides of the aquifer boundary.

The hydrogeochemical interpretation of the clusters further supports the spatial patterns of groundwater quality. Cluster 1, dominant in both inside and outside zones, represents baseline geogenic groundwater with relatively balanced major ions, indicating water–rock interaction under stable natural conditions and longer residence times (especially inside the aquifer). Cluster 2, which shows enrichment in Na⁺, Cl⁻, and TDS, reflects areas under anthropogenic influence such as irrigation return flows, or mixing zones with partially degraded water. It tends to occur in peripheral parts of the aquifer or near agricultural zones. Cluster 3, although less common, appears in both domains and may represent transitional waters, possibly originating from shared recharge sources or zones with similar lithological settings (e.g., carbonate-rich alluvial fans). The spatial distribution of clusters aligns well with geological heterogeneity and flow paths, indicating that HCA can effectively distinguish between natural and impacted water types.

PCA

The PCA reveals notable distinctions in groundwater chemistry between inside and outside the Neyshabur aquifer, reflecting differences in geochemical processes and aquifer conditions. The KMO values underscore better sampling adequacy inside the aquifer (0.79) compared to outside (0.44), indicating stronger correlations and suitability for multivariate analysis in the former (Table 5). The KMO value for the outside groundwater samples was slightly below the conventional threshold of 0.7, suggesting that the dataset may not be fully ideal for robust factor analysis. This indicates that the correlations among variables might not be strong enough to produce highly reliable factors. However, considering the exploratory nature of the analysis and the meaningful patterns identified, the factor analysis results were interpreted with due caution. Bartlett’s test is significant in both zones, validating the use of PCA. Inside the aquifer, three principal components were identified, accounting for 81.87% of the total variance with PC1 alone explaining 62.80%, PC2 at 9.75%, and PC3 at 9.30%. The high loading of major ions in PC1 (Ca²⁺: 0.86, Mg²⁺: 0.87, Na⁺: 0.79, Cl⁻: 0.91, SO₄²⁻: 0.95, EC: 0.93, TDS: 0.97, TH: 0.92) reflects the dominance of geogenic factors such as mineral dissolution and ion exchange processes within a chemically consistent aquifer. PC2 shows notable loadings for CO₃²⁻ (0.73) and pH (0.66) indicating the influence carbonate equilibrium and buffering capacity, particularly in areas with carbonate-rich formations. PC3, which includes loading on HCO₃⁻ (0.79) and pH (0.55), reflects variations in acid-base balance, potentially linked to CO₂ dissoltuion and recharge conditions. These components highlight local geochemical processes that are not solely driven by salinity but rather by carbonate weathering and aquifer-atmosphere interactions, especially in marginal and recharge zones.

Table 5. The Kaiser-Meyer-Olkin (KMO) and bartlett’s sphericity tests and the principal component analysis results for the study well and spring water samples

		Inside the aquifer			Outside the aquifer
KMO measure of sampling adequacy		0.79			0.44
Bartlett’s test of sphericity	Approx. chi-square	20,446			10,054
	df	66			66
	Sig	0			0
Rotation sums of squared loadings	Total	7.53	1.17	1.11	7.05	1.30
	% of variance	62.80	9.75	9.30	58.76	10.85
	Cumulative %	62.80	72.56	81.87	58.76	69.61
Number	Component	1	2	3	1	2
1	Ca2+	0.86	-0.05	0.19	0.77	0.34
2	Mg2+	0.87	-0.02	0.11	0.82	0.14
3	Na+	0.91	0.07	-0.10	0.94	-0.05
4	K+	0.76	0.04	-0.05	0.81	-0.07
5	CO32−	-0.21	0.73	-0.07	0.00	-0.72
6	HCO3−	-0.16	-0.35	0.79	0.04	0.64
7	SO42−	0.95	0.05	-0.06	0.94	-0.09
8	Cl−	0.90	0.02	0.01	0.84	0.20
9	pH	0.04	0.66	0.55	-0.03	-0.25
10	EC	0.97	0.03	-0.02	0.98	0.03
11	TDS	0.97	0.03	-0.01	0.98	0.03
12	TH	0.92	-0.05	0.18	0.79	0.29

Bold numbers refer to strong to moderate loadings

Extraction method: principal component analysis

Rotation method: varimax with Kaiser normalization

Outside the aquifer, only two principal components were extracted, explaining a lower cumulative variance of 69.61% (PC1: 58.76%, PC2: 10.85%), implying more heterogeneous hydrochemical conditions. While PC1 outside still shows strong loadings for Ca²⁺ (0.77), Mg²⁺ (0.82), Cl⁻ (0.90), SO₄²⁻ (0.84), EC (0.92), and TDS (0.93), the loadings are slightly reduced compared to the aquifer interior, indicating less uniform geochemical controls. Notably, HCO₃⁻ (0.64) appear in PC2 outside the aquifer, hinting at surface processes such as infiltration of HCO₃⁻-rich water. The absence of a third component outside the aquifer may reflect either the dominance of a few key processes or insufficient structure in the data due to more complex or variable recharge sources.

GWQI

The GWQI was separately calculated based on the results of the PCA and HCA analysis. For two domains, 9 physicochemical parameters with a loading weight of PC1 more than 0.7 (Ca²⁺, Mg²⁺, Na⁺, K⁺, SO₄²⁻, Cl⁻, EC, TDS, TH) were selected and the results are shown in Fig. 11. As shown in Fig. 11.a, the central and southeastern parts of the aquifer display large areas of very poor and poor groundwater quality, corroborated by the pie chart (Fig. 11.b), where 57.14% and 9.85% of the inside aquifer samples fall into these two classes, respectively. These degraded zones correspond to regions underlain by evaporative formations, such as gypsum and halite, which contribute to elevated concentrations of SO₄²⁻, Cl⁻, and TDS. In addition, prolonged over-extraction from deep wells has intensified water-rock interactions, leading to increased mineral dissolution and salinity. In contrast, groundwater quality outside the aquifer shows a relatively more favorable distribution, with 24.4% of samples categorized as good and 12.2% as excellent (Fig. 11.c). These areas are often located in the alluvial fan margins or higher elevation recharge zones, where lithology is dominated by coarse-grained, carbonate-rich sediments and there is more active infiltration from surface water and precipitation. The hydraulic connectivity between these zones and the main aquifer body appears limited, which may help preserve water quality outside the aquifer despite regional groundwater depletion. However, both zones still contain significant proportions of very poor water (32.2% outside the aquifer).

[See PDF for image]

Fig. 11

(a) Spatial distribution of different classes of GWQI; different classes of GWQI for the samples of inside the aquifer (b) and outside the aquifer (c)

Similar trends have been documented in arid and semi-arid settings globally. In the Nekor-Ghiss Plain of Morocco, high mineralization levels in groundwater were reported, primarily due to excessive withdrawals and drought-induced reduction in recharge (El Marnissi et al. 2023). In Sharjah, UAE, groundwater quality deteriorated markedly between 2004 and 2017, with salinity increases linked to seawater intrusion and limited natural recharge (Alhammadi et al. 2021). In contrast, studies in the Gurugram District of Haryana, India, showed better water quality in the central parts of the aquifer compared to its margins, where overuse and geological heterogeneity contributed to reduced quality (Kumar et al. 2024).

Conclusions

This study provides a detailed hydrogeochemical characterization of groundwater resources in the Neyshabur Plain, northeastern Iran, through an integrative approach combining classical graphical methods, multivariate statistical analysis, and spatial modeling. By distinguishing groundwater samples from within and outside the aquifer system, the research captures a comprehensive view of geochemical variability under contrasting hydrological and geological settings. The findings highlight that groundwater quality in the confined aquifer is predominantly governed by geogenic processes, including the dissolution of evaporitic minerals (e.g., halite and gypsum), prolonged water–rock interaction, and cation exchange, as evidenced by PCA and ion ratio analyses. In contrast, recharge zones outside the aquifer exhibit fresher water types dominated by carbonate dissolution, reflecting shallower flow paths and more active infiltration. Spatial interpolation of GWQI revealed alarming degradation in the central and southeastern parts of the aquifer, where over 70% of samples fall under poor to very poor quality classes, in stark contrast to relatively better conditions in external recharge zones. The hydrogeochemical classification, supported by Piper, Durov, and Gibbs diagrams, confirms a transition from Ca²⁺–HCO₃⁻ to Na⁺–Cl⁻ and Na⁺–SO₄²⁻ facies, tracking the chemical evolution of groundwater along flow paths influenced by lithological heterogeneity and evaporative concentration. The distinct clustering patterns derived from HCA and the strong factor loadings from PCA further underscore the spatial heterogeneity and complexity of groundwater quality dynamics across the plain.

Based on the findings, the following zone-specific management recommendations can be considered:

- Inside the aquifer:

1- Limit groundwater abstraction through regulatory controls and permit systems.

2- Prioritize artificial recharge projects in degraded zones, especially central and southeastern sectors.

3- Promote crop pattern changes and irrigation efficiency programs to reduce water demand.

- Outside the aquifer (recharge zones):

1- Legally protect recharge areas from contamination and uncontrolled land use.

2- Establish monitoring wells to detect early signs of quality deterioration.

3- Maintain low-impact agricultural practices to preserve groundwater quality.

To support decision-making, a three-tiered management framework is recommended:

1- Mapping-based zoning using GWQI to identify vulnerable areas.

2- Cluster-based prioritization, where groundwater types inform targeted interventions.

3- Policy integration, aligning water quality data with local planning, extraction licensing, and land-use regulation.

The integrated approach developed in this study can serve as a replicable model for groundwater monitoring and planning in other arid and semi-arid basins. Future studies could enhance this framework by incorporating socioeconomic data, land-use dynamics, and predictive modeling tools to support adaptive and sustainable groundwater governance.

Acknowledgements

I thank the Iranian Water Resources Management Company (IWRMC) for providing the valuable quality data of the well waters.

Author contributions

All the steps of manuscript preparation were done by the author.

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

Data availability

No datasets were generated or analysed during the current study.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

Abdalla FA, Ali RR, Fathy I (2019) Assessment of groundwater quality in the nile delta using multivariate analysis and GIS techniques. Environ Monit Assess 191(278). https://doi.org/10.1007/s10661-019-7427-3

Abdelaziz S, Gad MI, Tahan AH (2020) Groundwater quality index based on PCA: Wadi El-Natrun, Egypt. J Afr Earth Sci 172(103964). https://doi.org/10.1016/j.jafrearsci.2020.103964

Abtahi, M; Danesh, S; Davari, K; Qasemi, A. Investigation of Temporal changes in groundwater quality of Neyshabur plain and its possible causes. J Water Soil Conserv; 2015; 22, 4 pp. 171-186.

Ahmed, S; Rehman, MA; Khan, Z. Multivariate statistical techniques for the evaluation of groundwater quality in Pakistan. Hydrol Res; 2018; 2, 1 pp. 45-56. [DOI: https://dx.doi.org/10.1016/j.hydres.2018.05.001]

Alhammadi F, AlShamsi H, Abdelrahman M (2021) Groundwater quality degradation in sharjah, UAE (2004–2017): Temporal trends and geochemical indicators. E3S Web Conferences 234(01005). https://doi.org/10.1051/e3sconf/202123401005

Amiri Kaiz Kohani FS, Khashe Siuki A, Ahmadi Haji AH, Kavusi M (2021) Investigating the effect of groundwater on aquifer water salinity change (case study: Neyshabur Plain). Journal of Aquifer and Qanat 2(1):92–101. https://doi.org/JAAQ-1805-1007(R6)

Appelo, CAJ; Postma, D. Geochemistry, groundwater and pollution; 2005; 2 Amsterdam, A.A. Balkema:

Ariman S, Soydan-Oksai NG, Beden N, Ahmadzai H (2024) Assessment of groundwater quality through hydrogeochemistry using principal component analysis (PCA) and water quality index (WQI) in Kizilirmak delta, Turkey. Water 16(11). https://doi.org/10.3390/w16111570

Benyoussef S, Arabi M, El yousefi Y (2024) Assessment of groundwater quality using Hydrohemical process, GIS and multivariate statistical analysis at central rif, North Morocco. Environ Earth Sci 83(515). https://doi.org/10.1007/s12665-024-11798-6

Bouaissa M, Gharibi E, Ghalit M, Toupian JD, El Khattabi J (2021) Identifying the origin of groundwater salinization in the Bokoya Massif (central rif, Northern Morocco) using hydrogeochemical and isotopic tools. Groundw Sustain Dev 14(100646). https://doi.org/10.1016/j.gsd.2021.100646

Bourjila A, Dimane F, Ghalit M (2024) Appraising seawater intrusion in the Moroccan Ghiss-Nekor coastal aquifer: hydrochemical analysis coupled with GIS-based overlay approach. Desalination Water Treat 230(100612). https://doi.org/10.1016/j.dwt.2024.100612

Catania, ST; Reading, L. Hydrogeochemical evolution of the shallow and deep basaltic aquifers in tamborine mountain, Queensland (Australia). Hydrology J; 2023; 31, pp. 1083-1100. [DOI: https://dx.doi.org/10.1007/s10040-023-02617-6] 1:CAS:528:DC%2BB3sXlvFWrt78%3D

Chen W, Wu C, Pan S, Shi L (2024) Analysis on the Spatiotemporal evolutions of groundwater hydrogeochemistry and water quality caused by over-extraction and seawater intrusion in Eastern coastal, China. Front Earth Sci 12. https://doi.org/10.3389/feart.2024.1391235

Dao V, Urban W, Hazra SB (2020) Introducing the modification of Canadian water quality index. Groundw Sustainable Dev 11(100457). https://doi.org/10.1016/j.gsd.2020.100457

De la Hera Á, Rodríguez A, Sánchez J (2022) Hydrogeochemical characterization and groundwater classification in arid zones of Spain. Environ Earth Sci 81(4). https://doi.org/10.1007/s12665-0244-10123-4

El Marnissi B, El Faskaoui B, Maach A (2023) Hydrogeochemical assessment of groundwater in the Nekor-Ghiss plain, morocco: influence of drought and anthropogenic activities. Front Environ Sci 11(1179283). https://doi.org/10.3389/fenvs.2023.1179283

El Marnissi, B; El Faskaoui, B; Maach, A. Hydrogeochemical processes and groundwater quality in the central rif, Morocco. Environ Geochem Health; 2024; 46, 2 pp. 993-1009. [DOI: https://dx.doi.org/10.1007/s10653-023-01234-5]

El sheikh R, Soliman SA, Tawfik A (2021) Statistical evaluation of groundwater chemistry in the nile delta, Egypt. Groundw Sustainable Develoment 15(100656). https://doi.org/10.1016/j.gsd.2021.100656

Gad, A. Groundwater geochemistry and salinization processes in the Western nile delta, Egypt. Arab J Geosci; 2015; 8, 3 pp. 1501-1516. [DOI: https://dx.doi.org/10.1007/s12517-014-1260-9]

Gibbs, RJ. Mechanisms controlling world water chemistry. Science; 1970; 170, 3962 pp. 1088-1090.1:CAS:528:DyaE3MXjsVCmtw%3D%3D

Groeschke, M; Bosch, K; Daba, S; Kroemer, AL; Koeniger, P. Groundwater monitoring in challenging environments: an argument for the construction of observation wells based on data from niamey, Niger. Hydrology J; 2024; 32, pp. 1817-1831. [DOI: https://dx.doi.org/10.1007/s10040-024-02835-6]

Hosseinsarbazi A, Ismaeili K (2014a) Investigation and quantitative modeling of groundwater (case study: The Neyshabur Plain). Journal of Irrigation Science and Engineering 36(4):73–78. https://doi.org/20.1001.1.25885952.1392.36.4.8.5

Hosseinsarbazi, A; Ismaeili, K. Investigation of groundwater resources quality change on agriculture and technology (case study: the Neyshabur Plain). Iran J Irrig Drain; 2014; 8, 1 pp. 74-85.

Izeze EO, Imasuen OI, Badmus GO, Gbadebo AM (2023) Hydrogeochemical analysis and groundwater quality assessment of Ughelli south, Southern India. Environ Monit Assess 195(8). https://doi.org/10.1007/s10661-023-11580

Kanji S, Das S, Rajak C (2025) A comparative hydrochemical assessment of groundwater quality for drinking and irrigation purposes using different statistical and ML models in lower gangetic alluvial plain, Eastern India. Chemistry 372(144074). https://doi.org/10.1016/j.chemosphere.2025.144074

Khan, MA; Shakoor, MB; Nawaz, M. Groundwater quality evaluation using multivariate statistics in Southern punjab, Pakistan. Environ Sci Pollut Res; 2018; 25, 35 pp. 35754-35766. [DOI: https://dx.doi.org/10.1007/s11356-018-3381-2]

Khan MS, Farooq U, Aziz O (2019) Anthropogenic impact on groundwater chemistry in the Punjab region, Pakistan. Water 11(5). https://doi.org/10.3390/w11050890

Kumar, P; Singh, R; Meena, R. Groundwater quality evaluation in gurugram district, haryana, india: insights from geochemical and multivariate analyses. Environ Geochem Health; 2024; [DOI: https://dx.doi.org/10.1007/s10653-024-02047-8]

Li Y, Bian J, Li J, Ma Y, Anguiano JHH (2022) Hydrochemistry and stable isotope indication of natural mineral water in Changbai mountain, China. J Hydrology; Reg Stud 40(101047). https://doi.org/10.1016/J.EJRH.2022.101047

Mohammad, MAA; Szabo, NP; Mikita, V; Szucs, P. Taking the Spatiotemporal evolution of groundwater chemistry in the quaternary aquifer system of Debercen area, hungary: integration of classical and unsupervised learning methods. Environ Sci Pollut Res; 2025; 32, 11 pp. 6884-6903. [DOI: https://dx.doi.org/10.1007/s11356-025-36175-z] 1:CAS:528:DC%2BB2MXnvFWgsbs%3D

Mohseni U, Pande CB, Chandra Pal S, Alshehri F (2024) Prediction of weighted arithmetic water quality index for urban water quality using ensemble machine learning model. Chemosphere 352(141393). https://doi.org/10.1016/j.chemosphere.2024.141393

Naderianfar, M; Ansari, H; Ziaie, AN; Davary, K. Evaluating the groundwater level fluctuations under different Climatic conditions in the Neyshabur basin. J Irrig Water Eng; 2011; 1, 3 pp. 22-37.

Nagaraj, S; Masilamani, US. Hydrogeochemical and multivariate statistical approaches to investigate the characteristics of groundwater quality in fluoride-enriched hard rock region in Tirupathur district of Tamil nadu, India. Environ Sci Pollut Res; 2023; 30, 44 pp. 99809-99829. [DOI: https://dx.doi.org/10.1007/s11356-023-29254-6] 1:CAS:528:DC%2BB3sXhslKnt7bE

Ouarani M, Bahir M, Mulla DJ (2020) Groundwater quality characterization in an overall-located semi-arid coastal area using an integrated approach: case of the Essaouira basin, Morocco. Water 12(11). https://doi.org/10.3390/w12113202

Patel A, Rai SP, Akpataku KV (2023) Hydrogeochemical characterization of groundwater in the shallow aquifer system of middle Ganga basin, India. Groundw Sustain Dev 21(100934). https://doi.org/10.1016/j.gsd.2023.100934

Piper, AM. A graphic procedure in the geochemical interpretation of water analyses. Trans Am Geophysic Union; 1944; 25, 6 pp. 914-928.

Rahmati, S; Saadat, S; Poursoltani, MR. The effects of exploiting aquifers on the quality of groundwater resources (case study: Neyshabur plain, Northeast Iran). New Find Appl Geol; 2015; 17, 34 pp. 83-98.

Ramesh, R; Elango, L. Groundwater quality and its suitability for domestic and agricultural use in Tondiar river basin, Tamil nadu, India. Environ Monit Assess; 2012; 184, 6 pp. 3887-3899. [DOI: https://dx.doi.org/10.1007/s10661-011-2231-3] 1:CAS:528:DC%2BC38XmvVKjtLc%3D

Rao NS, Dinakar A, Sun L (2022) Estimation of groundwater pollution levels and specific ionic sources in the groundwater, using a comprehensive approach of geochemical ratios, pollution index of groundwater, unmix model and land use/land cover- a case study. J Contam Hydrol 248(103990). https://doi.org/10.1016/j.conhyd.2022.103990

Ravikumar, P; Somashekar, RK. Principal component analysis and hydrochemical facies characterization to evaluate groundwater quality in Varahi river basin, Karnataka state, India. Appl Water Sci; 2017; 7, pp. 745-755. [DOI: https://dx.doi.org/10.1007/s13201-015-0287-x] 1:CAS:528:DC%2BC2MXptFSitr4%3D

Shuaibu A, Kalin RM, Phoenix V, Lawal IM (2025) Geochemical evolution and mechanisms controlling groundwater chemistry in the transboundary Komadugu-Yobe basin, lake Chad region: an integrated approach of chemometric analysis and geochemical modeling. J Hydrology; Reg Stud 8(102098). https://doi.org/10.1016/j.ejrh.2024.102098

Sidibe AM, Lin X, Kone S (2019) Assessing groundwater mineralization process, quality, and isotopic recharge origin in the Sahel region in Africa. Water 11(4). https://doi.org/10.3390/w11040789

Singh SK, Kumar R, Tiwari A (2020) Hydrochemical evaluation of groundwater in semi-arid regions of Rajasthan. Appl Water Sci 10(5). https://doi.org/10.1007/s13201-020-01287-z

Taherian, P; Ansari, H; Davari, K; Beheshti, AA; Ziaei, AN. Modeling the groundwater quality (salinity) variations in Neyshabur plain using MODFLOW and MT3DMS. Iran J Irrig Drain; 2021; 15, 1 pp. 27-37.

Tahiri, M; Benaabidate, L; El Hmaidi, A. Hydrochemical characterization of groundwater using multivariate analysis: case study of the Souss-Massa basin, Morocco. J Afr Earth Sci; 2016; 121, pp. 75-87. [DOI: https://dx.doi.org/10.1016/j.jafrearsci.2016.04.011] 1:CAS:528:DC%2BC28XntlOqurw%3D

Tawfeeq, JMS; Disli, E; Hamed, MH. Hydrogeochmical evolution processes, groundwater quality, and non-carcinogenic risk assessment of nitrate-enriched groundwater to human health in different seasons in the Hawler (Erbil) and bnaslawa, Iraq. Environ Sci Pollut Res; 2024; 31, pp. 26182-26203. [DOI: https://dx.doi.org/10.1007/s11356-024-32715] 1:CAS:528:DC%2BB2cXmsFylsb8%3D

Umamageswari TSR, Thambavani DS, Liviu M (2019) Hydrogeochemical processes in the groundwater environment of Batlagundu block, Dindigul district, Tamil nadu: conventional graphical and multivariate statistical approach. Appl Water Sci 9(14). https://doi.org/10.1007/s13201-018-0890-8

World Health Organization (WHO). Guidelines for drinking water quality; 2017; Rome, WHO:

Yang, N; Su, C; Liu, W; Zhao, L. Occurrences and mechanisms of strontium-rich groundwater in Xinglong county, Northern china: insight from hydrogeological and hydrogeochemical evidence. Hydrology J; 2022; 30, pp. 2043-2057. [DOI: https://dx.doi.org/10.1007/s10040-022-02533-1] 1:CAS:528:DC%2BB38XisFWisr3P

Yazdani, V; Mansouriyan, H. The exploit potential zoning of groundwater resources by using of quantity and quality data in Neyshabur plain. J Irrig Water Eng; 2014; 4, 3 pp. 118-132.

Zeraati Neyshaburi, S; Pourreza Bilondi, M; Kashei-Siuki, A; Shahidi, A. Estimating the groundwater table of Neyshabur plain with introducing fuzzy possibilistic regression model. J Aquifer Qanat; 2022; 3, 1 pp. 53-64. [DOI: https://dx.doi.org/10.22077/jaaq.2018.1727.1008]