Combined Modeling of Multivariate Analysis and Geostatistics in Assessing Groundwater Irrigation Sustenance in the Middle Cheliff Plain (North Africa)

Abstract

The assessment of groundwater irrigation using robust tools is essential for the sustenance of the agro-environment in arid and semi-arid regions. This study presents a reliable method consisting of a combination of multivariate analysis and geostatistical modeling to assess groundwater irrigation resources in the Western Middle Cheliff (Algeria). For this goal, mean data from 87 wells collected during April to July 2017 were used. The hierarchical cluster analysis (HCA) using the Q-mode approach revealed three distinct water types, with mineralization increasing from cluster 1 to cluster 3. The Principal Component Analysis (PCA) utilizing the Varimax method approach allowed the extraction of three main components: the first and second (PC1, PC2), revealing that the geogenic process, have influenced the hydrogeochemical composition of groundwater. The pollution induced by agriculture activities has been related to PC3. Based on the combination of multivariate analysis and geostatistical modeling, the distribution maps were created by interpolating the factor distribution values acquired in the study region using the ordinary kriging (OK) interpolation method. The findings revealed that both natural processes and man-made activities have a substantial impact on the quality of groundwater irrigation. Cluster mapping, another often used combining approach, has shown its effectiveness in assisting groundwater resource management.

1. Introduction

During the last few decades, groundwater irrigation has gained a drastic amount of interest worldwide, especially in arid and semi-arid areas. On a global scale, almost 43% of the irrigated areas depend on groundwater resources [1]. In an agro-economy, the use of groundwater irrigation is moreover increasing in South Asia, the annual usage is 262 km³ yr⁻¹ (km³ = 1 billion m³), while it is (87 km³) in the Middle East, North Africa, and East Asia (57 km³). In terms of land cover, India (39 million hectares) has pioneered the use of groundwater in large areas, ahead of China (19 million hectares) and the United States (17 million hectares), which is in third place. In Algeria, the groundwater usage is 7.1 km³ yr⁻¹ and covered 88% of irrigated lands (883,004 ha) in the first decade of the twenty-first century [2]. Meanwhile, due to over-exploitation and the recent hardness of the climate, the groundwater quality becomes vulnerable to deterioration in these regions [3,4]. In addition, irrigation agriculture has proven to give better yield outcomes than rainfall agriculture within different environments. As a result, due to the highly competitive demand, efforts are needed to manage agricultural water as efficiently as possible [1,2].

Various pieces of research have assessed the sustainability of irrigated groundwater. In semi-arid regions, the geochemical properties on which it is based are controlled by several factors: anthropogenic activity, general geology, chemical weathering of several rocks’ varieties of water recharge quality, water-rock interactions, and ion exchange [5,6]. Therefore, there is a growing interest on the effect of geochemical processes on groundwater quality. However, considering the hydrogeochemical parameters, it is necessary to determine the suitability of groundwater for human needs and irrigation purposes [7]. Interestingly, groundwater quality related-studies in Algeria are now focusing on the irrigation prospects in various parts of the country [8,9,10,11,12,13,14,15,16,17,18,19].

The geochemistry of groundwater is often determined by its interaction with the mineral components of the flowing aquifer. Subsequently, many studies developed methods to characterize and assess the causes of the geochemical changes in groundwater. These methods have adopted approaches from multivariate statistical analysis and geostatistical techniques, or their combination [20]. Thus, the hierarchical cluster analysis (HCA) and principal component analysis (PCA) are very effective multivariate statistical techniques used to solve the issues coming from hydrological factors, mainly, the geochemical parameters [13,19,20,21,22,23,24], and geochemical composition controls [25,26,27] coming from the anthropogenic action [7,28,29,30]. The HCA, a multidimensional statistical technique, is a common tool in the classification and development of geochemical models based on observation data taken from factor values [31].

Recently, the use of Geographic Information Systems (GIS) has widened the horizon in determining the spatial distribution of groundwater irrigation quality and its mapping assessment using geostatistical techniques [13,19,32,33,34,35,36]. In Geostatistics, the kriging interpolation method allows one to predict data at non-sampled points and to map their spatial variability. Hence, Ordinary kriging (OK) is a very useful method in environmental studies, especially in groundwater quality related-studies [37,38].

The evaluation of groundwater irrigation quality is based on the physicochemical parameters compared to the international established limits [39] or some calculated indexes, i.e., sodium adsorption ratio (SAR). However, a successful management of groundwater irrigation requires the intervention of geostatistical tools using the kriging methods to gain better knowledge on the spatial distribution of the aquifer’s parameters. Nowadays, the available literature contains practical examples on the combination of kriging and factor analysis methods [7,40,41,42,43,44]. Moreover, the OK method is often utilized to interpolate factors representing the weights of specific processes that influence water chemistry. The study of these derived variables will reduce our multivariate workflow space and help to determine their spatial distribution in the submerged layer resulting from the calculation of variograms [7,40,41,42].

The Middle-Cheliff plain is one of the most important agricultural areas in Algeria, and indeed in North Africa, where there is a reliance on and very high demand for irrigation water. The intense agricultural activity comes from an imbalance between availability and the demand for groundwater. This study is based on an integrated investigation aiming to assess the groundwater quality in a semi-arid environment. The objectives of the study are as follows: (i) to identify the key hydrochemical features and major ion sources using a range of statistical analyses and hydrochemical techniques; and (ii) to identify the major ion sources using a variety of statistical analyses and hydrochemical methods; (iii) to assess the sustainability of groundwater for irrigation through the combination of multivariate analysis and geostatistical modeling.

2. Materials and Methods

2.1. Study Area

Towards the north-west of Algeria, the Western Middle-Cheliff Plain (WMCP) is located 200 km from the capital Algiers and 30 km from the Mediterranean Sea (Figure 1). The plain is surrounded by the Medjaja Mountains to the north, the Ourasenis Massif to the south, the Beni Rached Massif to the east, and the Lower Cheliff Plain to the west. It is located between 36°3′ and 36°15′ N and 01°4′ and 01°27′ E. It has a semi-arid climate with an average annual temperature of 19.4 °C, the summer temperatures are particularly high, reaching 29 °C in July. In winters, the temperature decreases significantly and it is about 10 °C in January (1980–2018). The study site receives an average of 320 mm from annual precipitation (Chlef City station).

The WMCP is an agricultural area, covering 27,000 ha, of 10,000 ha (37.4%) is reconquered for actual irrigation [45]. The plain is dominated by tree cropping (38.61%), especially orange trees; it is the most important irrigated cropping surface and it occupies the central peripheral part of the right bank of Oued Chlef. Vegetable cropping ranks in second place (34.56%) near the north borders (ABIADH-MEDJADJA, OULED-FARES, SOBHA, OUED-SLY, and OUM-DROU). Potatoes are the main crop, covering about 60% of the region. Cereals, with 26.56%, concentrated in the foothills, dominated by 90% of wheat and 10% of barley [45,46].

Surface water from the dams of Oued Fodda and Sidi Yakoub was historically the main source of water for irrigation systems in the investigated region. Nonetheless, long-term drought episodes have had an impact on irrigation canals. This led to the deterioration of the watering network. Thus, this was behind the exhaustion of the groundwater resource, and, more than before, it was necessary that this resource was reasonably managed. Given the importance of groundwater in irrigation practices [47,48], farmers from the study area have built well canals from groundwater to supply the cropping systems with new irrigation volumes. Therefore, groundwater resources became the main supply for irrigation in the studied plain starting from the early 2000’s [43].

2.2. Geological and Hydrogeological Contexts

The middle Cheliff watershed in Algeria corresponds to a subsided intramountainous furrow in the Tellian Atlas. It lies between the Dahra Massif in the north and the strongholds of the Ouarsenis in the south [49,50,51]. According to lithostratigraphic information [49,50], the Middle-Cheliff depression appeared as a whole during the Mio-Plio-Quaternary formation age. The Cheliff wadi, which enters the plain of Oum Drou (Pontéba) and exits at the Boukadir threshold, runs through them from east to west (Charon). These hydraulic thresholds correspond to the upwelling of the impermeable clay-marl substratum in the crossing, which effectively excludes any underground flow. The Middle-Cheliff plains’ intra-mountainous furrow was filled by Neogene deposits containing Quaternary, Pliocene, and Miocene sediments [49].

The Quaternary deposits are the starting point for the Neogene formations of marine origin, which can reach thicknesses of 3000 m [49]. They are most common in the plains, where they are composed of coarse alluvium (ancient Quaternary) and silt (late Quaternary) placed on upper Pliocene sandstone and limestone elements. The Lower Pliocene (Marine Pliocene) begins with a transgression on the Late Miocene gypsum series. Furthermore, it concludes with the Astian regression (Figure 2).

Three aquifers with varying hydrogeological potential exist in the region. The Upper Miocene limestone outcrops at the valley’s southern boundary and lies beneath the alluvium; the Pliocene sandstone, which is practically buried by the Quaternary formations; and the Pleistocene–Quaternary alluvial deposits, which constitute the valley’s embankment. Clays and marls coexist with sand, gravel, and conglomerate strata in these deposits. This last aquifer, which is the subject of our investigation, has an annual water withdrawal of around 15.5 hm³ [52].

The piezometric map created in 2012 [52] depicts hydro-isohypse closed curves in the center of the plain and open curves at the borders. This denotes a supply from the borders to the plain. Piezometric lines perpendicular to the valley axis intersect in the plain’s center before continuing in an East–West parallel route over the Oued Cheliff. The aquifer pumping does not compensate the weak recharge in the plain’s center.

2.3. Sampling Collecting and Analysis

Over 87 geolocated active wells were used for irrigation, and groundwater samples were collected throughout the study plain during the period from April to July 2017. During which, the irrigation of cropping systems extends all over the study area. The use of a GPS device helped in georeferencing the collected samples. To avoid sample contamination, the aquifers are pumped up for a few minutes before sampling. The samples were then put in preconditioned 1L polyethylene bottles and kept at 4 °C in an icebox. The pH and EC measurements were collected directly in the field with a multi-parameter handheld instrument WTW 350i (0% precision). The samples were then brought to the laboratory and examined for their major chemical contents (Ca²⁺, Mg²⁺, Na⁺, K⁺, Cl⁻, SO₄²⁻, HCO₃⁻, and NO₃⁻). The analytic methods utilized in this investigation matched the requirements of the American Public Health Association [53]. Calcium (Ca²⁺) and magnesium (Mg²⁺) concentrations were determined by volumetric measurements in the presence of EDTA aqueous solution; bicarbonate (HCO₃⁻) concentrations were determined by 0.1N hydrochloric acid titration; and chlorides (Cl⁻) concentrations were determined by titration of neutral medium solution of silver nitrates in the presence of potassium chromate. The amounts of sodium (Na⁺) and potassium (K⁺) were measured using flame photometry. In collaboration with a laboratory from the “Algerian Water Agency”, the spectrophotometry technique was used to doze sulfate (SO₄²⁻) and nitrate (NO₃⁻). The charge balance error (CBE) was assessed for each analyzed sample (Equation (1)). It was found to be acceptable in comparison to the (10%) range utilized in most laboratories [54].

$$\mathrm{Charge}-\mathrm{Balance} \mathrm{Error} \left(\mathrm{CBE}\right)=\left[\frac{{\displaystyle \sum \mathrm{cations}-{\displaystyle \sum \mathrm{anions}}}}{{\displaystyle \sum \mathrm{cations}+{\displaystyle \sum \mathrm{anions}}}}\right]\times 100$$

(1)

2.4. Statistical Analysis Approaches

The geochemical characterization of groundwater relies on the common use of the multivariate statistical analysis that retrieves valuable information from trends of chemical datasets [55]. In this study, the performance of HCA and PCA is usually used to extract the variability patterns and to identify the origins of the chemical elements. All the statistical processing was performed using the SPSS software v. 20.0 (IBM, Chicago, IL, USA) [56]. The HCA analysis used to classify the gathered groundwater samples is based on both Ward’s linkage approach and the Euclidean distance, both of which are frequently used in hydrochemical analyses [57]. The grouping process of samples with maximum similarity starts with the chemical compositions, then, a combination of the samples’ clusters with the Ward linkage method, the procedure is repeated until the full grouping of each sample. The clusters and the average linkage distance are presented graphically in the generated dendrogram from HCA to determine the number of homogeneous groups.

The PCA analysis creates associations between variables that reduces the sets of observations. The diagonalization of the correlation matrix is behind this reduction, giving rise to uncorrelated new dataset (orthogonal), with a decreasing order of importance called principal components (PCs) [58,59]. Considering the eigenvalue >1 [60], the execution of a Varimax rotation to these PCs has simplified the interpretation of the factors, whether it is related to the hydrochemical or anthropogenic processes that control the groundwater geochemistry. The factors with absolute loading value >0.75, 0.75–0.50, and 0.50–0.30 are considered “strong”, “moderate” and weak, respectively [61].

2.5. Geochemical Modelling

The geochemical modeling was performed by calculating the saturation index (SI) using PHREEQC v2.18 software for measuring data for defining the chemical reactions and aqueous speciation in the aquifer system [62]. The SI is calculated by Equation (2):

$${\gamma }_h=\frac{1}{2N(h)}{\displaystyle \sum _{i=1}^{N(h)}[Z(x_i})-Z(x_i+h){]}^2$$

(2)

where, IAP: ion activity product; K_sp: the solubility product at a given temperature.

The groundwater is being oversaturated when SI > 0 in a particular mineral (precipitation condition). SI < 0 indicates the under-saturation in a particular mineral (dissolution condition), whilst SI = 0 reflects the equilibrium state.

2.6. Geostatistical Modelling

Matheron was the first to develop and to apply the geostatistical theory starting from the mid-1960s [63]. Since then, it has grown into one of the most resilient approaches in applied statistics and earth sciences, particularly in hydrogeology, where it is used to create spatial maps of groundwater distribution and quality. The use of kriging interpolation techniques and semi-variogram models to interpolate the spatial pattern of groundwater irrigation analytical data is novel. The OK method, in particular, is an effective tool for linear interpolation [37,38].

The first step in the geostatistical modeling is determining the spatial dependence between neighboring observed points expressed using a variogram (γ_h) [64]. The semi-variance is defined as (γ_h), representing the semi-variance between the attribute values for all separated points by “h” distance as follows:

$${\gamma }_h=\frac{1}{2N(h)}{\displaystyle \sum _{i=1}^{N(h)}[Z(x_i})-Z(x_i+h){]}^2$$

(3)

where (γ_h): respresents the all pairs’ semi-variance for a lag distance h; Z(x_i): the groundwater irrigation parameter information per point i; Z(x_i + h): the groundwater irrigation parameter data for the rest of the separated points from x_i by a discrete distance h; x_i: the georeferenced location of Z(x_i); N(h): the number of observation pairs for separated points by a distance h [65].

The second step is the interpolation using the ordinary kriging method. The output consists of generated predictive maps representing the interpolation at unsampled locations of the regionalized variables with a minimum square error [66]. The calculation of the spatial distribution with ordinary kriging (OK) using the following equation:

$$Z^{\ast }(x_0)={\displaystyle \sum _{i=1}^n{\lambda }_i}Z(x_i)$$

(4)

and

$${\displaystyle \sum _{i=1}^n{\lambda }_i}=1$$

(5)

where, Z*(x₀) is the estimated value of the groundwater irrigation parameter by a location x₀, Z(x_i) is the available sample of the groundwater irrigation parameter per location x_i, and λ_i is the weight assigned to the sample value, and n is the number of considered samples in the prediction.

The combined use of multivariate statistical approaches and geostatistical modelling in this work aids in the evaluation of groundwater quality spatial mapping [20].

3. Results and Discussion

3.1. Hydrochemical Characterics of Groundwater

The concentrations for the majority of the chemical parameters were collected from nomalized data and were used to define the hydrochemical group using the Q-mod from the HCA approach, derived using Ward’s linkage method and the Euclidean distance for water sample similarity measurement. The HCA analysis produced a dendrogram (Figure 3), which divides groundwater samples into three groups. As a result, EC becomes a determining factor in identifying these clusters, ranging from cluster 1 to cluster 3. As a result, numerous clusters appear on the Piper [67] diagram to determine the geochemical development of groundwater types. In addition to statistical results, the physical and chemical properties of the three produced clusters were compared to Food and Agriculture Organization (FAO) criteria [39] and are given in

Table 1. Parameter’s value of the three principal water clusters.

Variables	Units	Cluster 1 (31 Samples)				Cluster 2 (47 Samples)				Cluster 3 (9 Samples)				Ayers and Westcot (1994) [39]
		Min	Max	Mean	SD	Min	Max	Mean	SD	Min.	Max.	Mean	SD
Ca	mg/L	25	198	109.6	38.9	73	292	196.4	49.9	186	500	366.2	32.8	400
Mg		4.74	89.5	38.6	21.1	40	241	111.6	45.9	67	415	263.4	33.1	60
Na		31.6	190	104.9	36.1	110	480	242.2	85.9	220	590	381.0	40.8	919
K		0.5	4.6	2.7	1.0	0	13	2.5	2.4	0	5	1.9	0.6	12
Cl		72	361	214.9	76.8	222.3	919.5	494.8	174.2	750	937.5	1129.8	113.7	1063
SO4		13	361	119.6	79.3	117.3	957	479.7	196.2	580	400	1077.8	97.6	960
HCO3		137	427	275.2	68.1	70	683.2	312.5	142.5	93	336	221.3	28.4	630
NO3		13.1	150	64.9	31.5	6	63.2	25.1	12.7	0	5	2.1	1.69	10
EC	µS/cm	802	1826	1358.7	282.9	1926	3658	2725.9	495.8	4550	8230	5843.3	1100.9	3000
pH	-	7.05	7.85	7.36	0.20	6.80	8.31	7.5	0.4	6.8	7.51	7.20	0.22	8.5

.

Cluster 1 relies on EC mean value equal to 1358.7 µS/cm measured from 31 wells, showing moderate water salinity for irrigation (C3) [68], this class can be used for cropping practices that are moderately salt tolerant. Given their abundance, the major ions order is: Ca²⁺ > Na⁺ > K⁺ and Cl⁻ > HCO₃⁻ > SO₄²⁻ (Figure 4). The hydrochemical type is distinguished by the evolution of Na⁺–K⁺–HCO₃⁻ into the Ca²⁺–Mg²⁺–Cl⁻–SO₄²⁻ facies (Figure 5). Bicarbonates (min = 137 mg/L, max = 427 mg/L, and mean = 275.2 mg/L) are the most abundant element in this cluster, followed by chloride (min = 72 mg/L, max = 361 mg/L, and mean = 214.9 mg/L) and SO₄²⁻ (min = 13 mg/L, max = 361 mg/L, and mean 119.6 mg/L). The majority of samples exceeded the acceptable nitrate level (10 mg/L) [39] for irrigation water (21 wells out of 31 exceeded the 10 mg/L irrigation guideline).

Cluster 2 is formed by 47 wells (54% of the total samples) characterized by high water salinity of 1947 < EC < 3673 µS/cm, mean = 2744.4 µS/cm. Sodium, chloride, and sulfates are the dominant ions. Meanwhile, the concentration of Na⁺, Cl⁻ and SO₄²⁻ varies from 110 to 480 mg/L, from 222.3 to 919.5 mg/L and 117.3 to 957 mg/L with mean concentrations of 242.2, 494.8 and 479.7 mg/L, respectively. The abundant major ions are, in descending order: Na⁺>Ca²⁺ >Mg²⁺ and Cl⁻ > SO₄²⁻- NO₃⁻ > HCO₃⁻ Figure 4). The hydrochemical type is characterized by Na⁺-K⁺-Cl⁻-SO₄²⁻ to evolve to Ca²⁺–Mg²⁺–Cl⁻–SO₄²⁻ facies (Figure 5).

In Cluster 3, only nine wells remain, with an EC value ranging from 4550 to 8230 µS/cm, mean = 5843.3 µS/cm, indicating highly mineralized groundwater. Mg²⁺ and Cl⁻ are the most dominant ions, stating a hyper composition of Ca²⁺–Mg²⁺-Cl⁻- SO₄²⁻ in the water facies. The influence of the environment on the wells of this cluster seems very important, especially from the agricultural practices that include chemical fertilizers and livestock.

For a better understanding of the hydrochemical process of groundwater, it is worth referring to Chadha (1999) [69], who presented a very clear illustration. This diagram is a modified version of the Piper diagram. It has been used successfully in various investigations to detect distinct hydrogeochemical processes [11,70]. The categorization of water using the Chadha diagram yielded four primary groupings (Figure 6). The application of this figure to groundwater data from the WMC plain reveals that the bulk of the three clusters’ samples fell into group 2, with fewer samples from cluster 2 falling into group 4. Group 2 is distinguished by an excess of Ca⁺² and Mg⁺² relative to Na⁺ and K⁺; this is due to the high calcium content of carbonate rocks, also the water-rock exchange along with the groundwater flowing (Figure 6). This group shows water type of Ca²⁺–Mg²⁺–Cl⁻–SO₄²⁻ composition.

3.2. The Principal Component Analysis (PCA)

The dataset of 87 observations that had undergone a PCA analysis was alongside ten (10) physicochemical parameters (EC, pH, Mg²⁺, Na⁺, Ca²⁺, K⁺, Cl⁻, SO₄²⁻, HCO₃⁻ and NO₃⁻). Through the SPSS software, the Varimax rotation method was performed (

Table 2. The score of PCA after Varimax rotation.

Variables	Component Data Values
	First Component	Second Component	Third Component
EC	0.955	0.006	0.093
Cl	0.911	−0.178	0.148
Ca	0.877	−0.023	0.212
SO4	0.867	−0.154	0.254
Mg	0.863	−0.205	0.290
Na	0.796	0.154	−0.055
pH	−0.060	0.983	−0.041
HCO3	−0.071	0.957	−0.195
K	−0.032	0.282	−0.837
NO3	0.394	0.049	0.769
Initial Eigenvalues of variances in %	48.087	20.861	15.578
Cumulative of variance %	48.087	68.948	84.526

). The results indicate complex geochemical components in groundwater. Following this, we extracted three factors with an eigenvalue >1, representing 84.526% of the total variance in the groundwater quality dataset (Table 2 and Figure 7).

According to Table 2, component 1 with the highest variance of 48.087% represents a strong mineralization of groundwater. Regarding, the positive important loading >0.75, the controlling parameters of the water table chemistry are mainly EC, Cl⁻, Ca²⁺, SO₄²⁻, Mg²⁺ and Na⁺ revealing salt leaching as well as a rock-water interaction within the aquifer. The major cations of this factor are Ca²⁺ and Mg²⁺, which may be formed by the weathering or reverse ion exchange process of carbonates [71].

The Rotated Varimax component 2 produced a total variance of 20.861% (Table 2), with a positive loading between pH and HCO₃⁻, with this component representing groundwater alkalinity. The pH in the examined groundwater is caused by the entry of precipitation into groundwater via deep percolation, as well as changes in mineral composition. HCO₃⁻ concentrations in groundwater revert to mineral dissolution in the saturated and unsaturated zones.

The Rotated Varimax component 3 represents 15.578% of the total sum of variation, with a significant loading of NO₃⁻ and K⁺, illustrating the influence of human-caused activities on the agriculture process. Indeed, the presence of nitrates in groundwater verifies the penetration of wastewater from agricultural practices that use a lot of nitrate fertilizers [11].

According to these findings, the key hydrogeochemical processes include aquifer material weathering, ion exchange, and human influences, notably agricultural activities.

The scatter-plot of both component 1 and 2 (Figure 7b) shows that all the water groups can be distinguished well in the PC space, and are the same as the clusters extracted from Q-mode HCA.

3.3. Hydrogeochemical Process

3.3.1. Evaporation

In semi-arid areas, the climatic effect through the evaporation process may increase the concentration level of all ions in the groundwater. The Na⁺/Cl⁻ vs. EC diagram plot is the indicator of evaporation (Figure 8a). The two chemical elements Na⁺ and Cl⁻ are moreover related to Halite dissolvance [72]. The ratio Na⁺/Cl⁻ < 1 indicates a Na⁺ content decrease related to its ion exchange with Ca²⁺ and Mg²⁺ in clays [73]. The increased ratio Na⁺/Cl⁻ (>1) refers to a none-halite source and indicates that the silicate weathering is the potential source of Na⁺. Figure 8a shows that the ratio of Na⁺/Cl⁻ lower than 1 (<1) is dominant and indicates that evaporation is the dominant process in the groundwater of the study area. Only four samples of cluster 1 present elevated values of Na⁺/Cl⁻ (12.9%), which indicates silicate weathering and 27 samples (87.1%) have a ratio >1 and the evaporation process. Fourteen samples of cluster 2 (29.9%) have a Na⁺/Cl⁻ ratio higher than one, and 33 samples (70.1%) have a Na⁺/Cl⁻ < 1. All the samples of cluster 3 have Na⁺/Cl⁻ ratio lower 1 and indicate that evaporation is dominated in this location.

The establishment of Plots between Na⁺ and Cl⁻ is helpful in controlling the salinization mechanism of groundwater. The halite dissolution gives rise to Na⁺ ions thus the resulting Na⁺/Cl⁻ molar ratio is close to one, indicating an obvious silicate weathering [74]. For the study area, the Na⁺ versus Cl⁻ plot (Figure 8b) is inferring that most of the points of clusters 1 and 2 are above the line of the freshwater evaporation and below the 1:1 section, and therefore are derived from a different anthropogenic background. The link between (Ca²⁺ + Mg²⁺) and (Na⁺ + K⁺) vs. the total cation concentration is typically used to confirm the silicate weathering process [11]. The association between calcium and magnesium ions and total cations reveals that fewer samples fall near to (Ca²⁺ + Mg²⁺) = 0.75 from the sum of cations line (Figure 8c), indicating that calcium and magnesium are much greater than potassium and sodium (Figure 8d). However, the majority of the calcium and magnesium in all clusters’ samples originates from the weathering of carbonate rocks.

3.3.2. Ion Exchange Process

To further study the processes affecting calcium concentration, (Ca²⁺ + Mg²⁺) versus HCO₃⁻ were used in groundwater samples (Figure 9a). For carbonates, the most common weathering reaction is simple dissolution, and the Ca²⁺/ HCO₃⁻ ratio is equal to 0.5. The low Ca²⁺/ HCO₃⁻ molar ratio of 0.5 (<0.5) indicates that the calcium and magnesium in the water are enriched by the exchange of sodium and/or the exchange of cations in the clay or the HCO₃⁻ that can come from the weathering of the silicate. In contrast, the high ratio (>0.5) suggests other sources of Ca²⁺ and Mg²⁺, such as reverse ion exchange, observed in hard rock formations with increasing salinity [75,76]. The water samples of the three clusters in the study area are located on the upper side of the 1:1 intersection, indicating that the excess of (Ca²⁺ + Mg²⁺) is behind the dominant carbonate weathering and contributing to the Ca²⁺ and Mg²⁺ groundwater [77,78]. Due to the breakdown of the quaternary calcareous layer, the molar ratio of (Ca²⁺ + Mg²⁺) ranges between 2.7 and 59.5. As a result of this, Ca²⁺ and Mg²⁺ must be balanced with ${\mathrm{SO}}_4^{-2}$ and Cl⁻ (Figure 9c).

The established plot (Ca²⁺ + Mg²⁺) versus (SO₄²⁻ + HCO₃⁻) defines three hydrogeochemical processes, mostly for samples at the 1:1 line, which are described by calcite and dolomite dissolution [79]. Figure 9b shows the distribution of the major samples on an equal line, indicating that the weathering degree of clay minerals, carbonates, and gypsum is relatively low, which is one of the mineralization factors of groundwater. The plot of Na⁺/Cl⁻ vs. Cl⁻ and (Ca²⁺ + Mg²⁺) vs. Cl⁻ (Figure 9d,e) shows that the salinity augments with the decrease in Na⁺/Cl⁻ and the increase in (Ca²⁺ + Mg²⁺), due to the reverse ion-exchange in the weathered clay layers [80]. Clay minerals have a boundary and negatively charged layered structure, on which the fixation and exchange of cations [81] occurs, as shown below:

$$\mathrm{Ca}\left(\mathrm{Mg}\right)-{\mathrm{Clay}}_2+2{\mathrm{Na}}^{+}\to {\mathrm{Ca}}^{2+}{(\mathrm{Mg}}^{2+})+\mathrm{Na}-{\mathrm{Clay}}_2$$

(6)

A decrease in the (HCO₃⁻)/(Cl⁻ + SO₄²⁻) ratio versus EC (Figure 9f) is related the elevation in electrical conductivity. The samples are grouped in two different poles in the study area: the first one describes the samples at their primitive state. It is dominated by HCO₃⁻ [11]. The second pole has an EC lower than 1500 µS/cm with the dominance of Cl⁻ and SO₄²⁻ ions characterizing an anthropogenic activity. The samples of this pole have EC > 1500 µS/cm, where human activities contribute to mineralization from many sources (agricultural practices, discharge of untreated sewage, etc.). This pole is dominated by samples from Cluster 2 and 3.

The creation of the bivariate plot from [(Ca²⁺ + Mg²⁺) − (HCO₃⁻+ SO₄²⁻)] versus [(Na⁺ + K⁺) − Cl⁻] helped in studying the relationship between the ion exchange and reverse ion exchange processes [82]. In the present study, the obtained results are showing a controlled action of the reversed ion exchange process on the chemistry of the groundwater with a slope value of −1 (Figure 9g).

3.4. Geochemical Modelling

The determination of the groundwater geochemistry is based on the interacted action between the water and the aquifer’s mineral components. The use of SI to predict the reactive rock mineralogy from groundwater sampled data, eliminating the need to collect solid samples and to avoid mineralogy analysis. In this study,

Table 3. Summary of minerals Saturation Index (SI) in groundwater using PHREEQC.

	Anhydrite	Calcite	Dolomite	Gypsum	Halite
Cluster 1
Min	−3.11	−1.02	− 0.87	−1.57	−2.89
Max	−1.21	0.91	1.05	2.05	−0.99
Mean	−1.84	0.08	0.22	0.29	−1.62
SD	0.39	0.39	0.39	0.70	0.39
Cluster 2
Min	−1.76	−0.89	−0.75	−1.37	−1.54
Max	−0.77	1.63	1.77	3.6	−0.55
Mean	−1.14	0.39	0.53	1.15	−0.92
SD	0.23	0.64	0.64	1.23	0.23
Cluster 3
Min	−1.10	−0.72	−0.58	−0.92	−0.88
Max	−0.53	0.54	0.68	1.42	−0.31
Mean	−0.71	0.08	0.23	0.63	−0.50
SD	0.17	0.42	0.43	0.80	0.17

presents the outputs of the SI calculation for the Calcite, Aragonite, Dolomite, Gypsum, Anhydrite, and Halite selected minerals. Positive SI values indicate oversaturation and mineral to precipitation from groundwater. The SI negative values refer to an under-saturation and thus mineral dissolution in the groundwater. The state of equilibrium occurs when the saturation index is between −0.5 and +0.5 and indicates that the mineral does not dissolve or precipitate in this groundwater [83]. The SI results indicate a precipitation tendency of carbonate minerals (calcite, aragonite and dolomite) in all clusters (Figure 10a). Given the high evaporation and low precipitation rates (<400 mm/year) in the study area, this is a reasonable record characterizing a deposition of aragonite, calcite and dolomite [13]. However, anhydrite, gypsum and halite (Figure 11b,c) are in a state of under-saturation, which states that the concentrations of their soluble components Na⁺, Cl⁻, Ca²⁺ and SO₄²⁻ do not require a mineral equilibrium [84], while anhydrite and gypsum minerals are on the way to equilibrium. Due to the nature of Na⁺, it exists in low concentrations compared to Cl⁻, which can, on one hand, associate with clay minerals through an ion exchange process. On the other hand, the decrease in concentration of Na and increase in concentration of Ca returns mainly to the reverse ion-exchange that led to the minimized dissolution of gypsum [85].

3.5. Assessment of Groundwater Irrigation Parameters

The adequacy of groundwater for irrigation is determined by comparing the risks of salinity and sodium concentration to the total cations in the system. The sodium concentration is expressed as a percentage of sodium (% Na) and a ratio of sodium adsorption (SAR), whilst salinity is expressed as electrical conductivity (EC) (

Table 4. Sustenance of groundwater for irrigation based on EC, %Na and SAR.

Parameters	Range	Class	Number of Samples	Percent
EC (µS/cm) (Wilox 1955) [68]	C1 < 250	Excellent	0	0
	250 < C2 < 750	Good	0	0
	750 < C3 < 2250	permissible	41	47.13
	2250 < C4 < 5000	Doubtful	39	44.83
	C5 > 5000	Unsuitable	7	8.04
% Na (Wilox 1955) [68]	0–20	Excellent	5	5.75
	20–40	Good	63	72.41
	40–60	Permissible	15	17.24
	60–80	Doubtful	4	4.6
	>80	Unsuitable	0	0
SAR (Richards 1954) [87]	S1 < 10	Excellent	87	100
	10 < S2 < 18	good	0	0
	18 < S3 < 26	Doubtful	0	0
	S4 > 26	unsuitable	0	0

All cations are expressed in mmol_c/L.

). High salinity water tends to decrease the osmotic potential of irrigated plants over time. Meanwhile, the high concentrations in Na⁺ ion participates in a cation exchange process that affects the soil’s physical properties [86]. Many studies have shown that irrigation water enriched with Na⁺ leads to cation exchange reactions with a release of magnesium and calcium into the water. On the other hand, the use of water enriched with Ca²⁺ for irrigation leads to the release of Na⁺ ions and the preferential adsorption of Ca²⁺ ions at the cation exchange sites.

To determine the aptitude of groundwater for irrigation, EC is a very important key-parameter (Table 4). Hence, around 47.13% of the samples were admissible for irrigation in the study area (class C3). More than fifty samples were ranged from Doubtful class (C4: 44.83%) to Unsuitable class (C5: 8.04%), and no samples were Excellent class or good class (C1 and C2) for irrigation.This result shows that the high EC value affects the productivity of the plants (deficit in transpiration) and the structure of the soil. As a result, irrigation water with high EC rates diminishes crop production potential [68]. The sodium percentage (%Na) indicated that 72.41% of the samples are good quality for irrigation, while 17.24% are permissible for irrigation and 4.6% have doubtful quality. Meanwhile, five samples (5.75%) were considered as excellent quality. SAR values are in the range below 10, which according to SAR indicates excellent water and could be used for most soil types.

The irrigation parameters of %Na an SAR are expressed by following equations:

$$\%\mathrm{Na} =\left(\frac{\mathrm{Na}+\mathrm{K}}{\mathrm{Ca}+\mathrm{Mg}+\mathrm{Na}+\mathrm{K}}\right)\times 100$$

(7)

$$\mathrm{SAR}=\frac{\mathrm{Na}}{\sqrt{\frac{(\mathrm{Ca}+\mathrm{Mg})}{2}}}$$

(8)

Usually, the evaluation of water quality for irrigation is given using both a US salinity diagram and a Willox diagram (Figure 11); the first, the US salinity diagram, was proposed by the USSL [87], to show that the sodium adsorption ratio (SAR) has been plotted against the EC (Figure 11a). It is established using the DIGRAMME software.

According to Figure 11a, all the groundwater samples of cluster 1 belong to the categories C3S1. It is average to poor quality used with precautions and requires drainage, leaching and/or gypsum application. The samples of cluster 2 were ranged C4S1 (61.7%), C4S2 (34%), and C4S3 (4.3%). The quality of most samples (C4S1) is poor to bad water quality, to be used with caution for heavy soils and sensitive plants; the samples of C4S2 is qualified as very poor-quality water, and was used only for light and well-drained soils and for resistant plants with the necessity of leaching doses and/or gypsum contribution; the samples of C4S3 are very poor quality for irrigation. Most of the samples of cluster 3 (66.7%) fall into the field of C5S2; not recommended for irrigation.

The Willox diagram [68], which connects salt percentage and electrical conductivity values, demonstrates that, with the exception of one sample, the bulk of the groundwater samples are of excellent quality. Cluster 2 samples were 72.3% full in the dubious to inadequate (26.7 %) quality range. All of the samples in Cluster 3 are of poor quality (Figure 11b). The high sodium content in irrigation water causes it to be absorbed by clay particles, displacing Mg²⁺ and Ca²⁺ ions. The exchange of Na⁺ in water for Ca²⁺ and Mg²⁺ in soil affects soil permeability [88].

3.6. Geostatistical Modelling

The irrigation parameters of the studied water have been used for performing the Kriging interpolation. The best-fitted semi-variogram models were selected based on the nugget’s variance/sill ratio. There were three forms of spatial dependence considered in this study: strong when the ratio is less than 25%, mild when the ratio is between 25% and 75%, and weak when the ratio is greater than 75%. The mean error (ME) and root mean square standardized error (RMSSE) data are used to calculate the precision of the forecasts. The ME should be near zero, and the RMSSE should be close to one [38,89]. To assess the correctness of the interpolation, the cross-validation method was used. During cross-validation, each measured point is gradually deleted, and the value is predicted using the remaining data. The difference between each measured and predicted value is the error. Cross-validation can also be used to determine the best option [90].

According to

Table 5. Best-fitted semi-variogram models and cross-validation for OK of groundwater irrigation parameters and factors analysis.

	Variogram Parameters						Prediction Precision
Parameters	Transformation	Model	Nugget (C0)	Partial Sill (C)	$(\frac{{\mathbf{C}}_0}{{\mathbf{C}}_0+\mathbf{C}})\times 100 (\%)$	Spatial Dependence	Mean	Root-Mean- Square Standardized Error
EC	Log	Exponential	0.0154	0.2927	4.99	Strong	0.125	0.9441
%Na	Log	Rational Quadratic	0.0499	0.0471	51.44	Moderate	−0.178	1.0283
SAR	Log	Spherical	0.0783	0.1506	34.21	Moderate	−0.009	1.0440
Factor 1	Original Data	Spherical	0.0968	0.9916	8.89	Strong	0.0132	1.1242
Factor 2	Original Data	Spherical	0.8051	0.2252	78.14	Weak	−0.0073	0.9972
Factor 3	Original Data	Spherical	0.1388	0.8247	14.40	Strong	−0.0362	1.0304

, the EC parameter was fitted to the exponential semi-variogram model, %Na fitted Rational to the Quadratic model and SAR were feted to the Spherical model. Depending on the nuggets’ variance/sill ratio, EC has a strong structure of spatial dependence, while %Na and SAR have a moderate spatial structure. All three parameters were interpolated after transforming the dataset on logarithmic value.

The spatial distribution of the EC parameter (Figure 12a) shows the high salinity risk by groundwater irrigation. The water of doubtful quality (C4) is dominant, in particular in the south and northeastern part of the study area. Although, the north east of the study area has less EC values belonging to the permissible class (C3). The C5 class, is located in the center and a little to the west of the area and relies on very high salinity risk (EC > 5 dS/m). This class (C5) is located when there is excessive groundwater pumping (Figure 2), which, together with evaporation, is most likely the reason. The overlay analysis for %Na spatial distribution (Figure 12b) indicated that the upper part of the north-west and lower south-west zones had a higher concentration of %Na. This result is similar to the SAR distribution in the study area (Figure 12c). The lower concentration of SAR was observed mainly in the north-east part of the study area at the localities of Ouled Fares and Labidh Medjadja. Furthermore, the groundwater indicated calcium and magnesium precipitation in the form of calcite and dolomite (SI > 0). This can result in a reduction in the molality of Ca²⁺ and Mg²⁺. Due to evaporation, the value of SAR rises [91].

Combination of Factor Analysis and Geostatistical Modeling

In the current study, the estimation of spatial variability of the dominant factors retrieved their scores in the study area. The factor scores utilized in the study were taken from the PCA factor score coefficient matrix using standardized data for each water sample site. The OK technique is used to build spatial variability maps by an interpolation of the factor scores for each sampled point [20,44]. The spatial variability maps are generated. All the factor analysis was fitted to the Spherical model (Table 5). Factors 1 and 3 have strong spatial distribution dependence with nuggets’ variance/sill ratio <0.25. Factor 2 has weak spatial distribution dependence; the ratio is more than 0.75.

Figure 13a displays the distribution of Factor 1 scores in the study area. The central portion of the study area represents the high positive scores. The EC distribution map and the Factor 1 distribution pattern are genuinely similar for the study area (Figure 13a). This confirms that the interaction of the evaporation-water process in the aquifer is the determining factor of the overall groundwater chemistry in the study area and is behind the high scores of Factor 1.

The map representing the scores of factor 2 (Figure 13b) indicates high values to the northwestern and towards southeastern parts of the study area where the geologic processes enrich the aquifer with more carbonate and, thus, higher pH values in this location.

The factor 3 distribution map (Figure 13c) is under the effect of natural and anthropogenic processes with an increase in values from the west to the east in the study area. All the agricultural activities are mostly practiced in the center and eastern part of the study area (Figure 13c).

Another representation of groundwater quality from combining HCA and spatial analysis is cluster mapping using Mapping Clusters tools of ArcGIS software. This representation can show the locations of the cluster in the study area. Figure 13d shows the spatial distribution of the three clusters in the area. The groundwater quality for irrigation increasing from cluster 1 to cluster 3, the spatial distribution shows that cluster 2 is dominant in the study area, the cluster 1, qualified by good quality for irrigation, is located in the north area at the localities of Ouled fares and the southwest (BouKadir). Cluster 3, high risk, is located in the center and a small part in the north of Boukadir district to the extreme west of the WMCP.

Finally, these results have serious consequences on the environmental and agriculture sustenance of the study plain. Considering the increasing potential use of groundwater resources for irrigation, it will certainly lead to an important degradation risk of soil quality and endangering the crop yield, both in the short and long terms.

4. Conclusions

The triple combination of hydrogeochemical process, multivariate analysis and geostatistical modeling proved its usefulness in assessing the sustainability of irrigation groundwater in the agricultural plain of Western Middle-Cheliff, Northwest Algeria.

Three water types were extracted using a Q-mode cluster analysis based on the groundwater quality data sets. The first cluster represents the low salinity (EC = 1358.7 µS/cm) water sample and stretches to the north of the study area. The second one (Cluster 2) is dominant in the area. Also, a domination of Ca–Mg-Cl- SO₄²⁻ facies in this cluster indicates a questionable salinity class C4 (EC = 2744.4 µ/cm). The third group (Cluster 3) has a high salinity level, is marked as unsuitable for irrigation (EC = 5843.3 µ/cm) and is dominated by Ca-Mg-Cl and SO₄²⁻.

The anionic concentrations of Cl⁻ and SO₄²⁻ dominant in groundwater, as well as Ca²⁺ and Mg²⁺ from cationic abundances, show that alkaline earth (Ca²⁺ and Mg²⁺) outnumber alkali metals (Na⁺ and K⁺), and strong acids (Cl⁻ and SO₄²⁻) outnumber weak acids (HCO₃, CO₃⁻). The hydrogeochemical mechanisms that influence groundwater chemistry have been dominated by reverse ion exchange. For evaporating minerals, all water groups are under-saturated. Carbonate minerals, on the other hand, are over-rated in all clusters.

The Varimax and PCA approaches produced three PCs that are said to account for more than 84 percent of the total variance. Both PC1 and PC2 discovered that the geogenic process has attacked the groundwater hydrogeochemical composition, in which the dissolution of carbonate and evaporated rocks, as well as reverse ion exchange and weathering processes, are included. The PC3 gene is connected to farming operations, in which irrigation becomes more frequent during the plant’s development phase.

Using ordinary kriging of geostatistical analysis, the geographical distribution maps for groundwater irrigation parameters indicated a significant spatial structure for EC and a moderate spatial structure for Na and SAR. The EC distribution map reveals that the uncertain quality (C4) predominates, particularly in the research area’s south and northeast. The class C5 is located in the plain’s center and is impacted by the region’s heavy pumpage as well as the evaporation process. The precipitation of Ca²⁺ and Mg²⁺ can raise the value of SAR (>10) in the investigated region.

The combination of PCA and the geostatistical modeling to assess groundwater suitability for irrigation was based on the factor scores in each groundwater monitored well. They are injected as variables in the OK method creating various surfaces and displaying the range and the rate of groundwater quality influenced by these common factors. The spatial distribution of Factor 1 seems strong and it is similar to the EC parameter spatial distribution in the area. This EC factor characterizes the groundwater mineralization. Factor 3 distribution reflects the recourse to irrigation at the actual time of the plant’s development. The obtained cluster map shows their spatial distribution along the area.

The pairing of geostatistical modeling and multivariate statistical methods in the evaluation of groundwater sustenance for irrigation provides a powerful and effective tool for studying and analyzing common factors of uncertainty and completing any long-term monitoring data analysis, particularly for datasets collected at large-scale regions. The current findings will aid future studies on environmental clean-up, pollution prevention, natural variation, and the implementation of water management initiatives.