Introducing a new geostatistical approach to classify groundwater samples based on Stiff diagram: Case study of Chahardoly aquifer, west of Iran

Abstract

Groundwater quality is pivotal for sustainable resource management, necessitating comprehensive investigation to safeguard this critical resource. This study introduces a novel methodology that integrates stiff diagrams, geostatistical analysis, and geometric computation to delineate the extent of a confined aquifer within the Chahrdoly aquifer, located west of Hamadan, Iran. For the first time, this approach combines these tools to map the boundaries of a confined aquifer based on hydrochemical characteristics. Stiff diagrams were used to calculate geometric parameters from groundwater chemistry data, followed by simulation using a linear model incorporating the semivariogram parameter γ(h). The Root Mean Square Error (RMSE) of the linear model was used to differentiate confined from unconfined aquifers based on hydrochemical signatures. Validation was conducted by generating a cross-sectional hydrogeological layer from well logs, confirming the presence of aquitard layers. The results successufully delineated the confined aquifer's extent, showing strong agreement with hydrogeological log data. By integrating stiff diagrams with semivariogram analysis, this study enhances the understanding of hydrochemical processes, offering a robust framework for groundwater resource identification and management.

Graphical abstract

Keywords

Geostatistics / Stiff diagram / Semivariogram / Confined aquifer / Chahardoly / Asadabad

Introduction

Groundwater reservoirs are vital water resources that support various human activities, including agriculture, industry, and drinking water supply. Their significance is particularly pronounced in semi-arid and arid regions, where surface water is scarce. The hydrochemical characteristics of groundwater are crucial for assessing its suitability for domestic, agricultural, and industrial use (Chegbeleh et al. 2020). Numerous factors influence groundwater quality, including natural processes such as rock-water interactions and weathering (Shuaibu et al. 2024; Aweda et al. 2023; Eyankware and Omo-Irabor, 2019), as well as anthropogenic activities. Scientists have employed various approaches to study hydrochemistry, including geostatistical and statistical techniques (Afzal et al. 2024), to understand the interaction between geological formations and water quality (Yidana et al. 2012).

Previous studies have extensively applied statistical methods to analyze groundwater chemistry. Principal Component Analysis (PCA) has been widely used to reduce large hydrochemical datasets and identify dominant factors affecting groundwater composition (Wold et al. 1987; Krishan et al. 2023). Correspondence analysis, another statistical method, is often utilized to measure the similarity between water samples and categorize hydrochemical types (Ji et al. 1995). Factor Analysis (FA) has also been employed to investigate key hydrochemical components, revealing that mineral dissolution and silicate weathering are major contributors to groundwater quality variations (Sun, 2014).

Geostatistical methods such as semivariogram analysis and kriging techniques have been applied to understand the spatial distribution of groundwater quality parameters (Krasowska and Banaszuk, 2017; Farzaneh et al. 2022). Ordinary Kriging (OK) is a widely accepted interpolation method for estimating groundwater characteristics in unmeasured locations (Masoud, 2014). Studies incorporating semivariogram analysis have demonstrated its effectiveness in mapping hydrochemical variations and delineating contamination zones (Sa´nchez-Martos et al. 2001).

Hydrochemical diagrams, including Piper, Gibbs, and Durov diagrams, have been extensively used to classify groundwater types based on ion concentrations. Stiff diagrams, in particular, offer a visual representation of groundwater composition, yet their quantitative potential remains underexplored (Sun, 2014; StiffJr, 1951). While Stiff diagrams are commonly applied for qualitative hydrochemical assessments, few studies have attempted to integrate them into a quantitative framework.

The Chahardoly confined aquifer, located in western Iran, plays a critical role in supporting agricultural and domestic water demands. Understanding its hydrochemical behavior is essential for sustainable groundwater management. This study introduces an innovative approach by integrating Stiff diagram geometric properties with semivariogram analysis to quantitatively assess groundwater hydrochemistry. Unlike traditional applications of Stiff diagrams, this method enhances groundwater classification by incorporating spatial correlation analysis. By applying this novel framework, we aim to improve the accuracy of confined aquifer delineation and contribute to more effective groundwater resource management.

1 Materials and method

1.1 Area of Chahardoly Plain

A plain covering an area of 89 km² is located in the west of the Hamadan Province. The Chahardoly Plain lies between latitudes 48°02'N and longitudes 34°50′E. The Shahab River, which flows from the surrounding mountains into the Asadabad Plain, is the most important river in the area. East of the Chahardoy Plain lies the confined aquifer, covering around 6.5 km² and adjacent to Chenar village (Moradi Nazar Poor et al. 2022). Limestone and metamorphosed rock formations are located to the east of the plain, alongside the confined aquifer. The dry season spans from June to September, characterized by minimal precipitation, while the wet seasons extends from October to May, with the highest rainfall (Fig. 1 and Fig. 2). The average annual precipitation over the past 44 years has been 336 mm.

1.2 Explanation of the Chahardoly confined aquifer

The Chahardoly confined aquifer is bounded by two aquitard layers that separate it from the underlying aquifer. The features of the unconfined aquifer are evident, such as three artesian wells and tow springs. Examples of such natural water sources include artesian wells and springs (Moradi Nazarpoor and Jafari, 2019).

1.3 Methodology

This study proposed several methodological frameworks to clarify processes involved. The steps have been summarized as follows:

(1) Calculating geometric characteristicsof all Stiff diagrams;

(2) Determining the boundary area of the confined aquifer;

(3) Validating the results.

The methodology is schematically presented step by step in Fig. 3.

1.3.1 Groundwater Sampling

To collect data for analysis, 18 groundwater samples from both the confined and unconfined Chahardoly aquifers were collected with appropriate distribution, ensuring sufficient representation in both aquifers. The samples were then transferred to Shahrood University's water laboratory. Using titration, pH and EC meters, and a chromatography machine, pH, Electrical Conductivity (EC), and the major ions (Calcium, Magnesium, Sodium, Potassium, Chloride, Sulfate, and Bicarbonate) were measured (Fig. 4).

To verify the accuracy of the major ion measurements, the charge balance error between cation and anion concentrations were calculated (Yidana et al. 2012) using Eq. 1. The charge balance error formula (Appelo and Postma, 2004) was applied to ensure that the total cations and anions were within 10% of each other (Table 1).

(1)$ \text{Error}\text=\frac{\displaystyle\sum \text{cation}-\displaystyle\sum \text{anion}}{\displaystyle\sum \text{cation}+\displaystyle\sum \text{anion}}\times 100 $

1.3.2 Step1 (Geometric characteristics)

This step comprises two main parts. First, one sample is designated as the primary sample, while the remaining samples were compared to it. In this study, W7, located within the confined aquifer, is selected as the primary sample. As an artesian well (Fig. 4), W7 represents the confined aquifer. Second, the geometric properties of Stiff diagrams, such as line lengths, are determined by the intervals between ions on the Stiff diagram.

To calculate the interval between nodes, as illustrated in Fig. 5, it is necessary to scale the ion values dimensionless, treating them as distance distances. Subsequently, the geometric characteristics of each sample are compared with those of the primary sample. During this comparison, the Root Mean Square Error (RMSE) between the primary sample and each other sample is calculated. The sample with the lowest RMSE is considered the most similar to the primary sample (Eq. 2).

(2)$ RMSE=\left[\frac{1}{n}\sum \left(\widehat{x}-{x}_{i}\right)\right] $

1.3.3 Step2 (Boundary area determination)

To estimate the boundary of the confined aquifer and differentiate samples belonging to the confined aquifer from those of the unconfined aquifer, a semivariogram was utilized to quantify the dissimilarity between the two groups of samples. If it is assumed that samples closer in proximity are more similar–a relationship that can be analyzed using the semivariogram (Johnston et al. 2001). The covariance structures of measured sample points with respect to distance is characterized by the semivariogram (Moonchai and Chutsagulprom, 2020).

The experimental semivariogram, γ(h), is calculated as follows (Eq. 3).

(3)$ \gamma \left(h\right)=\frac{1}{2NP\left(H\right)}\sum _{i=1}^{NP\left(h\right)}\text{‍}\left[Y\right({x}_{i}+h)-Y({x}_{i}){]}^{2} $

The number of pairs of samples separated by h is denoted as N(h). With i and x_i represent a sample point. At sampling site x_i, γ(x_i) is the observed value of the variable (Delhomme, 1978). Several wells exhibiting the highest similarity to W7, representative of the confined aquifer, were selected based on their depth and RMSE values (Table 2). Subsequently, the optimal liner model was determined using the γ(h) values of these wells (Fig. 6), and all experimental semivariogram values were used to calculate the theoretical value (Table 3).

An assumption was made that in homogeneous areas (i.e. the confined aquifer is considered a homogeneous area), a linear model exists where the value of γ(h) increases steadily as the distance between samples increases. The homogeneous model line was determined and drawn after plotting the semivariogram (Figs. 6 and 7).

To determine the optimal boundary between the two aquifers, both the theoretical γ(h) (estimated by the homogeneous line) and the experimental γ(h) values (calculated from sample data), were compared. The error between these values was calculated, and based on well logs analysis, samples with an error less than 25% were classified as part of the confined aquifer.

(4)$ Error=\frac{\gamma (h{)}_{t}-\gamma (h{)}_{e}}{\gamma (h{)}_{t}}\mathrm{*}100 $

1.3.4 Step3 (Validation)

To validate the proposed method, the results from a previous study on the Chahardoly, based on well logs, were reviewed (Moradi Nazarpoor and Jafari, 2019). Rockwork software was used to generate a cross-section of the geological layers for this case study. The extent of the unconfined aquifer was determined by identifying the aquitard layer through geological layer simulation. Finally, the boundary of the unconfined aquifer obtained from the proposed method was compared with the results derived from geological layer simulation.

2 Results and discussion

2.1 Result of the geometric characteristic

Six parameters for anions and cations are required to construct a stiff diagram for each sample. Table 1 presents the chemical analysis results of the Chahardoly Aquifer, encompassing 18 samples collected from both confined and unconfined aquifers, along with the concentrations of major anions and cations. Fig. 4 shows the locations of these samples.

For each sample (X1 to X13), Table 2 provides the geometric characterization results (i.e. the interval between nodes) in 13 sections, along with the Root Mean Square Error (RMSE) calculated between the primary sample and others. Notably, the RMSE values for W16, W18, W30, and W26 are higher than those of other samples. This indicates that as spatial distance increases, the RMSE tends to increase as well. As expected, the RMSE between W7 and itself is zero due to perfect similarity.

2.2 Results of the variogram

The variogram of the samples was analyzed to investigate the relationship between the confined and unconfined aquifers. Both the experimental values (i.e. measured variogram values for each sample) and theoretical values (i.e. variogram values calculated based on linear model assumption) of γ(h) are illustrated in Fig. 6 and Table 3.

Fig. 7 displays the value of (γ(h)) plotted against the sample distances. The value of the linear model alongside the theoretical and experimental γ(h) calculations are summarized in Table 3.

The linear model, representing the homogeneity of the confined aquifer, assumes a linear increase in γ(h) values with increasing distance. This indicates that as the distance between samples increases, the variogram value γ(h) also increases linearly.

In summary, the linear model was employed to calculate the theoretical values for each sample, while experimental values were derived from actual observations. Samples exhibiting low error values demonstrate a close alignment with the linear model, indicating strong similarity between experimental and theoretical values. For instance, the primary sample, W7 shows similar positioning to W27, W14, and W13. Conversely, samples W16, W30, W18, and W12 display significantly higher error values, corresponding to their greater distance from W7. Based on error value for W3, it is likely located outside the confined aquifer, within the unconfined aquifer. Additionally, the error values for W9, W10, and W11 show a decreasing trend, while the errors for W5 and W6 remain relatively high, indicating that these also belong to the unconfined aquifer (Fig. 8). Furthermore, all aquifer points, whose values were calculated using the simplest interpolation method (Spline), were plotted alongside the samples in Fig. 8.

2.3 Validation based on well logs

Fig. 9 shows the boundary of the confined aquifer, which lies east of the Chahardoly aquifer. Three hydro-geology cross-sections are presented, namely A-A, and B-B. These cross-sections are based on the results of the study by Moradi Nazarpoor and Jafari (2019). Fig. 9 illustrates the condition of the aquitard layers and the confined aquifer.

Depending on the elevation of the aquitard layer and the locations of artesian wells and springs, a confined aquifer can be clearly identified in a specific area. The confined aquifer covers approximately 6.5 km². Aquifer simulations and well logs were used to establish the hydrogeological layers (Moradi Nazarpoor and Jafari, 2019). Using their results, an error range was assessed to determine the maximum allowable error, which was set at 25% for identifying wells within the confined aquifer.

Fig. 10 shows that the extent of the confined aquifer was determined using both the proposed method and well log data. The size and location of the aquifer are similar for both approaches. The yellow squares in Fig. 10a represent locations identified as unconfined aquifer by the proposed method, while Fig. 10b displays the unconfined aquifer area based on hydrogeological well logs.

On the others hand, the confined aquifer samples (W7, W24, W14, W27, W13, and WS1) generally exhibit relatively lower sulfate (SO₄²⁻) and chloride (Cl⁻) concentrations compared to the unconfined aquifer samples, However, their hydrochemical profiles appear similar, likely due to the predominantly carbonate geological formation. This suggests a reduced influence from surface contamination and anthropogenic sources such as agricultural runoff. Bicarbonate (HCO₃⁻) concentrations are slightly higher in the confined aquifer, indicating a greater influence of carbonate weathering and longer water-rock interaction times. For example, W7 (2.6 Meq/L) and W14 (2.3 Meq/L) show relatively high HCO₃⁻ values. Magnesium (Mg²⁺) and calcium (Ca²⁺) levels are moderate in both aquifers, but confined samples generally show a more stable composition, reflecting limited external influences. Sodium and potassium (Na⁺ + K⁺) concentrations are relatively consistent across both aquifers, with minor variations. However, confined samples show slightly lower variability, which may indicate less impact from evaporation and surface water mixing (Table 4).

2.4 Assumptions and Limitations of the new proposed method

While this study relies on certain assumptions, such as the adequacy of the sampling network to clearly define the primary sample as representative of the confined aquifer and other water sources for comparison, the accuracy and precision of ionic measurements are also critical. Ensuring a well-distributed and comprehensive sampling network is essential for accurately capturing the hydrochemical characteristics of the aquifer.

3 Conclusion

This study introduces a novel methodology that integrates geometric analysis of Stiff diagrams with semivariogram analysis to delineate confined aquifer boundaries in the Chahardoly aquifer, Iran. By extracting geometric characteristics from Stiff diagrams across thirteen intervals, the approach effectively characterized major ions in groundwater samples. Semivariogram analysis distinguished between confined and unconfined aquifer samples, while a linear model, with an error criterion of less than 25%, defined the confined aquifer boundary. Validation through hydrogeological cross-sections confirmed that samples with errors below 25% aligned with the confined aquifer, closely matching well log data. This innovative method enhances groundwater quality assessment, crucial for supporting domestic, agricultural, and industrial needs in semi-arid regions. However, recognizing the limitations of a single approach, future research could integrate geophysical surveys or machine learning to further refine aquifer boundary delineation.

References

[1] Afzal M, Liu T, Bao AM, et al. 2024. Groundwater hydrodynamics and supply sustainability in Punjab, Pakistan: A geological statistical approach using the GAMEAS algorithm based on data from 1105 boreholes. China Geology, 7: 1−16. DOI: 10.31035/cg2024066.

[2] Appelo CAJ, Postma D. 2004. Geochemistry, groundwater and pollution. CRC press. DOI: 10.1201/9781439833544

[3] Aweda AK, Jatau BS, Goki NG. 2023. Groundwater geochemistry and hydrochemical processes in the egbako aquifer, northern bida basin, nigeria. RBRH, 28: e31. DOI: 10.1590/2318-0331.282320230010.

[4] Chegbeleh LP, Akurugu BA, Yidana SM. 2020. Assessment of groundwater quality in the talensi district, northern ghana. The scientific world Journal, 2020(1): 8450860. DOI: 10.1155/2020/8450860.

[5] Delhomme JP. 1978. Kriging in the hydrosciences. Advances in Water Resources, 1(5): 251−266. DOI: 10.1016/0309-1708(78)90039-8.

[6] Eyankware M, Omo-Irabor O. 2019. An integrated approach to groundwater quality assessment in deter- mining factors that influence the geochemistry and origin of sandstone aquifers southern niger delta region of nigeria. Malaysian Journal of Geoscience, 23–32. DOI: 10.26480/mjg.02.2019.23.32.

[7] Farzaneh G, Khorasani N, Ghodousi J, et al. 2022. Application of geostatistical models to identify spatial distribution of groundwater quality parameters. Environmental Science and Pollution Research, 29(24): 36512−36532. DOI: 10.1007/s11356-022-18639-8.

[8] Ji H, Zhu Y, Wu X. 1995. Correspondence cluster analysis and its application in exploration geochemistry. Journal of Geochemical Exploration, 55(1-3): 137−144. DOI: 10.1016/0375-6742(95)00025-9.

[9] Johnston K, Ver Hoef JM, Krivoruchko K, et al. 2001. Using ArcGIS geostatistical analyst, volume 380. Esri Redlands.

[10] Krasowska M, Banaszuk P. 2017. Geostatistical methods in the assessment of the spatial variability of the quality of river water. In E3S Web of Conferences, 22: 00089. EDP Sciences. DOI: 10.1051/e3sconf/20172200089.

[11] Krishan G, Bhagwat A, Sejwal P, et al. 2023. Assessment of groundwater salinity using principal component analysis (PCA): A case study from Mewat (Nuh), Haryana, India. Environmental Monitoring and Assessment, 195(1): 37. DOI: 10.1007/s10661-022-10555-1.

[12] Masoud AA. 2014. Groundwater quality assessment of the shallow aquifers west of the nile delta (egypt) using multivariate statistical and geostatistical techniques. Journal of African Earth Sciences, 95: 123−137. DOI: 10.1016/j.jafrearsci.2014.03.006.

[13] Moonchai S, Chutsagulprom N. 2020. Semiparametric semivariogram modeling with a scaling criterion for node spacing: A case study of solar radiation distribution in thailand. Mathematics, 8(12): 2173. DOI: 10.3390/math8122173.

[14] Moradi Nazar Poor S, Jafari H, Safari M, et al. 2022. Hydro-geochemistry of the asadabad plain, west of Iran: Application of statistical methods. Geopersia, 12(2): 369−378. DOI: 10.22059/geope.2022.339940.648653.

[15] Moradi Nazarpoor S, Jafari H. 2019. Determining the confined aquifer in Chahardolly (Asadabad). 3rd Iranian National Conference on Hydrology, University of Tabriz.

[16] Sa´nchez-Martos F, Jimenez-Espinosa R, Pulido-Bosch A. 2001. Mapping groundwater quality variables using pca and geostatistics: A case study of Bajo Andarax, Southeastern Spain. Hydrological Sciences Journal, 46(2): 227−242.

[17] Shuaibu A, Kalin RM, Phoenix V, et al. 2024. Hydrogeochemistry and water quality index for groundwater sustainability in the Komadugu-yobe Basin, sahel region. Water, 16(4): 601.

[18] Stiff Jr HA. 1951. The interpretation of chemical water analysis by means of patterns. Journal of Petroleum Technology, 3(10): 15−3.

[19] Sun L. 2014. Statistical analysis of hydrochemistry of groundwater and its implications for water source identification: A case study. Arabian Journal of Geosciences, 7: 3417–3425.

[20] Wold S, Esbensen K, Geladi P. 1987. Principal component analysis. Chemometrics and Intelligent Laboratory Systems, 2(1-3): 37−52.

[21] Yidana SM, Ophori D, Banoeng-Yakubo B, et al. 2012. A factor model to explain the hydrochemistry and causes of fluoride enrichment in groundwater from the middle voltaian sedimentary aquifers in the northern region, ghana. ARPN Journal of Engineering and Applied Sciences, 7(1): 50–68.

RIGHTS & PERMISSIONS

Journal of Groundwater Science and Engineering Editorial Office