Groundwater quality assessment using revised classical diagrams and compositional data analysis (CoDa): Case study of Wadi Ranyah, Saudi Arabia

2.1 Location and climatic context

The governorate of Ranyah is part of the Mecca region, located in the western part of Saudi Arabia, with a total area of around 62,000 km². It lies between latitudes 42° 30′ and 43° 00′ E and longitudes 21° 00′ and 21° 30′ N (Fig. 1). This region is characterized by a semi-arid climate, with cold winters and hot summers, with clear seasonal variations in rainfall (average 450 mm/year) and high evaporation (average 400 mm/year) (Subyani and Al Ahmadi, 2010; Saleem et al., 2020). According to the study realized by Geng and Boufadel (2017), which examined the effects of evaporation and precipitation rates on groundwater salinity and dynamics using field data and numerical models, indicates that groundwater salinity in the surface layer increases with increasing evaporation rate. Another important point is that not only climatic factors affect water resources, but also human and mining activities such as fertilizers, wastewater, agricultural effluents and industrial waste, all of which generate anthropogenic sources of pollution, while rock weathering also generates geogenic sources of pollution.

Close

Fig. 1

Geological map and inventory of groundwater samples.

Export to PPT

According to the preliminary results of the general population and housing census carried out by the Ministry of Statistics and Information (2018), 47,034 people live in the Ranyah region. The majority of the governorate's inhabitants are employed in agriculture, with the Ranyah Valley offering excellent conditions for this activity. The valley is considered an important water reservoir supplying the region, and lies within the Hedjaz Mountains, part of the Arab Shield that runs north–south parallel to the Red Sea. Wadi Ranyah rises in the Baha Mountains, where its altitude varies from 1,750 to 1,880 m above mean sea level (amsl), and extends eastwards into the Saudi Arabian desert, where its altitude gradually decreases (Mohamed et al., 2023).

2.2 Geology setting

The project plan to map the geography and geology of the Arabian Peninsula received approval in July 1953 from the governments of the Kingdom of Saudi Arabia, the United States Department of State, and the Arabian American Oil Company (ARAMCO), with neighboring countries also granting their approval. The main goal was to make the fundamental geology of the area, as mapped by ARAMCO and the U.S. Geological Survey, accessible to the public. This collaboration between governments, oil companies, and private individuals marked a significant milestone in geological research (Brown et al., 1989).

The specific study area is in the western part of KSA, on the Arabian Shield. This shield covers about a third of the Arabian Peninsula, or about 630,000 km² of land (32 % of the Kingdom's total area). It consists predominantly of Precambrian rocks, mainly igneous and metamorphic rocks, with younger volcanic rocks found in the western region (Brown et al., 1989). Fig. 1 provides a visual representation of the various rock formations from the Precambrian to the Quaternary that make up the geology of Wadi Ranyah. These rocks include andesite and basalt flows, granodiorite, diorite, granitic plutons, pyroclastic dacite and rhyolite, schist, gneiss, and metamorphosed volcanic wacke and sandstone. The area's Quaternary alluvial layers, on the other hand, consist of sands, silts, clayey sands, silts and gravels that were eroded from nearby soil and collected in Wadi Ranyah. Their thickness ranges from 5–10 m in the upper parts to 10–15 m or more in the lower parts (Mohamed et al., 2023).

Hydrogeologically, the Wadi Ranyah aquifer consists of the alluvial deposits of the wadi, with a thickness ranging from 3 m upstream to 12 m downstream, as well as in the fractured bedrock. The depth of the water table fluctuates between 3 and 10 m, with no systematic variability along the course of the wadi (Subyani and Al Ahmadi, 2010).

3.1 Sample collection and analysis

In this study, sampling focused on shallow wells to extract groundwater from Quaternary alluvium and fractured rock in the Ranyah region. 45 samples were collected from urban and rural areas after 25–30 min of pumping time (Fig. 1). High-density polyethylene bottles were washed with hydrochloric acid and rinsed 2 or 3 times with sample water before sampling. Each bottle was meticulously numbered and labelled before being transported. Samples collected after being filtered with quantitative Whatman Grade 44 filter paper and acidified with nitric acid (pH < 2) were stored in a cooler at 4 °C before being transferred to the laboratory for analysis of physicochemical parameters (APHA, 2017). The pH and TDS of the site were all measured using a pH meter (HANNA – HI 991300) which was calibrated before reading the results. In addition, each groundwater sample was analyzed for major ions (SO₄⁻², HCO₃^–, Cl⁻, NO₃^–, Ca⁺², Mg⁺², Na⁺, K⁺, Cl⁻) and minor ions (F, Mn, Pb and Zn) in milligrams per Liter. Titrimetric and photometric methods with the HANNA spectrophotometer (HI 83215) were employed to analyze the main elements: sulfates, bicarbonates, chlorides, magnesium, calcium, sodium and potassium. While trace elements of groundwater samples were determined using inductively coupled plasma mass spectrometry (ICP-MS).

All measurements and analysis procedures mentioned above were carried out in accordance with the WHO international standard. The analytical error (less than ±10 %) of the analyzed ion content was verified by the charge balance error (CBE) percentage based on the following Eq. (1):

3.2 Compositional data analysis (CoDA)

Many researchers, particularly in geosciences, are interested in geochemical data as a subset of multivariate data, which belongs to the category of compositional data and consists of a matrix of nonnegative relative values with a constant sum (Egozcue et al., 2003; Greenacre, 2018; Reghais et al., 2024). In detail, CoDA, which is based on obtaining relative information by comparing the compositional components, was developed by Aitchison (1982), and Egozcue et al. (2003). According to Egozcue et al. (2003), and Greenacre (2018), the use of traditional statistical techniques can lead to false assumptions and correlations. To address this problem, two general approaches were later developed.

• The first approach, known as “ staying in the simplex”, with a constant sum K and a non-negativity constraint (Equation (2), consists of a constraint on a set of D components or part within a D-1-dimensional hyperplane in positive real space (R+D − 1), called “simplex space” (Egozcue et al., 2003).

• The log-ratio transformations are used in the second approach to transfer the simplex into real space, where conventional Euclidean geometry and standard statistical methods can be used. This is achieved using either the centered log-ratio transformation (clr) or the isometric log-ratio transformation (ilr), as in the works of Aitchison (1982), and Egozcue et al. (2003).

3.3 Revised diagrams (Compositional perspective)

Numerous graphical and numerical methods have been developed in the literature to elucidate the intricate nature of geochemical characterization of groundwater. This complexity arises from a multitude of uncertainties stemming from natural factors, such as geological structures, aquifer reservoirs, and transit time, as well as anthropogenic factors like wastewater discharges and agricultural activities (Reghais et al., 2023). In general, conventional analysis methods such as ion ratios, correlation analysis, and graphs like the Gibbs diagram have been used to achieve this goal. However, these techniques often overlooked the relative structure of data sets, leading to inaccurate and incorrect results (Wheeler et al., 2021). Consequently, it is now crucial to update schematics, while still enabling straightforward interpretation of the geochemical findings (Templ et al., 2008).

Recently, new statistical tools have been developed to effectively compare and handle digital information. One of these techniques involves converting data into ilr coordinates, which allows for easy transfer from simplex space to real space. This technique has been proven to be valuable both in theory and practice, especially in the construction of balances using the Sequential Binary Partition (SBP) of composition x. The SBP divides composition parts into two non-overlapping groups at each of the D-1 partition a step, resulting in ilr coordinates of D-1. Equation (5) represents the equilibrium between these groups in RD-1 (Pawlowsky-Glahn and Buccianti, 2011; Buccianti, 2015).

In the ith step of SBP, x_i+ and x_i- represent the groups of parts that are separated. r_i+ and r_i- indicate the number of parts in ci+ and ci−, respectively. Additionally, the function g(.) corresponds to the geometric mean of its argument.

This paper utilizes two classical diagrams that have been recently revised in light of an approach called compositional data analysis. These diagrams are used to describe the groundwater in Wadi Ranyah, Saudi Arabia. The first diagram is the ilr-plot, which is a new type of diagram that incorporates the information found in Piper's diagram. The ilr-plot was introduced by (Shelton et al., 2018). The second diagram is a revised version of the Gibbs diagram, which takes into account the properties of compositional data. This revised version was introduced by (Buccianti, 2015). These diagrams allow for the exploration and description of the chemistry of water in an exploratory and descriptive manner. The balances of the modified and replaced diagrams are presented for each axis (Supplementary Materials Table S2), as well as the sequential binary partition for the ilr diagram (Supplementary Materials Table S3).

3.4 Drinking water resources assessment

Groundwater plays a crucial role in Saudi Arabia's overall water supply. However, recent studies by Alghamdi et al. (2023) indicate that the quality of these resources has been deteriorating due to excessive groundwater extraction in the dry region. This degradation poses significant challenges to human health. In this section, our focus lies on analyzing chemical data from 45 groundwater samples, encompassing 15 chemical parameters including Ca²⁺, Mg²⁺, Na⁺, K⁺, Cl⁻, HCO₃^–, SO₄²⁻, NO₃^–, TDS, EC, F, Pb, Mn, and Zn. To comprehensively assess the hydrogeochemical condition of groundwater sources with regards to their potability, we utilized the weights derived from the relative loadings and corresponding eigenvalues of the water parameters in the principal components (PCs). These weights were employed in the calculation of the Water Quality Index (WQI) as per Equation (6), as described by Ewuzie et al. (2021).

The weight of parameter i, noted Wi, was obtained from the eigenvalue of each PC as well as the factor loadings of each parameter based on the PCA results. Ci and SI correspond to the value of the parameter measured for each sample and the WHO (2017) guideline, respectively. The water quality index (WQI) remains among the reliable methods for assessing water quality according to drinking water guidelines (Shanmugamoorthy et al., 2023). It categorizes drinking water quality into five classes: Excellent quality when WQI ≤ 50, Good quality (51 ≤ WQI ≤ 100), Poor quality (101 ≤ WQI ≤ 200), Very poor quality (201 ≤ WQI ≤ 300) and, finally, Unsuitable for consumption when WQI > 300 (Meng et al., 2016).

4.1 Multivariate compositional statistics (CoDA-MSA)

4.1.3

4.1.3 Robust principal component analysis

Log-transformed ratios are important for assessing compositional data because they emphasize the relative magnitudes and variability of the components rather than their absolute values (Aitchison, 1982). To ensure proper statistical analysis of the variable inputs obtained from the laboratory, a log-transformation (clr) was applied prior to conducting the principal component analysis. The dataset consists of 45 groundwater samples and 15 variables, including 8 major ions (Mg²⁺, Ca²⁺, K⁺, Na⁺, SO₄⁻², Cl⁻, HCO₃^–, and NO₃^–) and 4 minor ions (F, Zn, K, Pb, and Mn), in addition to TDS, EC, and pH. However, the variation to be explained determines the number of robust components that can correctly represent the data. Consequently, the two robust PCs were selected (rPC 1: 37.8 %, and rPC 2: 20.1 %), representing a total variance of 57.9 %. The descriptive statistics obtained from the principal component analysis of the clr-transformed data are presented in Supplementary Materials (Table S5). Based on the quartile and median values, it is evident that the second cluster exhibits lower PC1 values compared to the first cluster (I). Nevertheless, it is worth noting that cluster I still contains samples with high PC1 values (>6). Conversely, cluster II consists of a sub-cluster with high PC2 values (>4.5).

The F1-F2 factorial plan (57.9 % of variance) reveals an opposition between two groups of variables (Fig. 2). The F1 axis, explaining 37.8 % of the variance, represents the main differentiating factor between these groups. The first group, referred to as G1, consists of pH, manganese, and trace elements. These variables may indicate the weathering effect of crystalline rocks and agricultural activities. Notably, this group corresponds to wells located in the upstream part, specifically the Al-Hujrah zone, which is considered the recharge zone. On the other hand, the second group (G2) can be subdivided into two subgroups: the subgroup (G2-1) corresponds to the wells from the Ranyah zone, which are rich in HCO₃^–, Cl⁻, K⁺, SO₄²⁻ and Na⁺. While the subgroup (G2-2) represents the downstream wells in the Al-Hujrah zone, which are mineralized like those in Ranyah, but richer in calcium and nitrates.

Close

Fig. 2

Principal component analysis of clr-transformed data: Biplot of the PC1 versus PC2 factors.

Export to PPT

Regarding axis F2 (20.1 % of variance), it distinguishes mineralized waters rich in HCO₃^–, Cl⁻, K⁺, SO₄²⁻ and Na⁺ that are well correlated with electrical conductivity (G2-1), which corresponds to the wells in the Ranyah zone, from highly diluted waters (G1) with higher concentrations of trace metallic elements.

4.1.4

4.1.4 Cluster analysis

Hierarchical cluster analysis (HCA) in R mode was applied to the 15 hydrochemical parameters, while HCA in Q mode was utilized to examine the 45 samples of groundwater. The resulting dendrograms enabled the classification of the waters in three distinct groups (G1, G2, and G3) (Fig. 3). This sequential logic (G1-G2-G3) aligns with the flow direction, as it corresponds to the concentration of elements along the pathways and flow lines of groundwater.

Close

Fig. 3

Dendrograms: (A) cluster analysis of variables; (B) cluster analysis of samples.

Export to PPT

The first group, referred to as G1, consists of trace metal elements, pH, TDS, and EC. These characteristics are used to define wells in the recharge zone, where high contents of lead, zinc, and fluoride are present. These elevated levels are a result of the weathering of igneous and metamorphic rocks, as well as agricultural activities (Subyani and Al Ahmadi, 2010).

The second group, known as G2, focuses on calcium (Ca⁺²) and nitrate (NO₃^–) levels (Fig. 3A). This group includes five wells: H51, H55, H58, H61, and H62, which are located downstream of the Al-Hujrah zone. The water from these wells is mineralized and contains high levels of calcium. The excess calcium and nitrates compared to other wells is likely due to agricultural activities.

The third group, G3, includes all wells in the Ranyah area, which is the discharge zone. This group exhibits a strong correlation among major ions, namely Mg⁺², HCO₃^–, Na⁺, SO₄⁻², K⁺, and Cl⁻ (Fig. 3B). The waters in this group are mineralized, with an electrical conductivity (EC) greater than 750 μS/cm. This mineralization is attributed to the alteration of silicate minerals, such as plagioclase, pyroxenes, and amphiboles, which are commonly found in igneous rocks like andesite, granodiorite, and mafic igneous rocks (Drouiche et al., 2022). Additionally, the evaporation effect and overexploitation of the groundwater aquifer for irrigation (Reghais et al., 2023).

4.2 Isometric log-ratio-ion plot

Although the Piper diagram has been a useful tool for categorizing different types of water for the past eight decades, it has also presented challenges. The diamond-shaped structure of the Piper diagram renders it difficult to differentiate between fluids primarily comprising of Ca⁺², Mg⁺², Cl⁻, and SO₄⁻², especially when dealing with large data sets. To address this issue, a new version of the Piper diagram has been developed. This updated diagram incorporates compositional data analysis (CoDA) to generate four (4) isometric log-ratio (ilr) coordinates for the six (6) principal chemical parameters analyzed in the traditional Piper diagram. These coordinates are then utilized to create a four-field biplot, referred to as the ilr-ion plot. The ilr-ion plot allows for the visualization of additional data, such as the relationships between Ca²⁺/Mg²⁺ and Cl⁻/SO₄⁻². It effectively conveys information from the traditional Piper diagram regarding composition, water chemistry types, and geochemical processes (Shelton et al., 2018).

The information in Piper's diamond can be understood by referring to the anion triangle (SO₄⁻² + Cl⁻) versus (HCO₃^–) and the cation triangle (alkaline-earth vs. alkali), as shown in Fig. 4. Based on the top right panel of the diagram (Fig. 4B), bicarbonates (HCO₃^–) are generally the dominant anions in all three groups of water samples, with an average content of about 170 mg/L. This aligns with the HCO₃^– content found in the crystalline domain. A slight increase in HCO₃^– content was observed downstream, which can be attributed to variations in rainfall. The river experiences dry periods of more than three months, followed by cycles of heavy rainfall, leading to precipitation and a higher level of HCO₃^– in the groundwater (Saleem et al., 2020).

Close

Fig. 4

Isometric log-ratio-ion plot for Ranyah groundwater samples.

Export to PPT

All samples from the third group (Ranyah zone) exhibited negative values, indicating a higher abundance of alkaline elements (Na⁺ + K⁺) compared to alkaline-earth elements (Ca⁺² + Mg⁺²) in relationship to the flow direction (Fig. 4C). In the top right panel, data points for groups (1) and (2) have negative y-coordinate values and positive x-coordinate values (Fig. 4). The sampling points within the cation panel represent the Ca⁺² type (Fig. 4C), while in the anion panel, they represent the HCO₃^– type (Fig. 4B), suggesting an HCO₃^– − Ca⁺² type of water. The relative abundances of Ca⁺² and Mg⁺² (x-axis) compared to SO₄⁻² and Cl⁻ (y-axis) can be distinguished in the left top panel (Fig. 4A). In this panel, all samples from the three groups have positive x-coordinates, indicating that the Ca⁺² content was mostly higher than the Mg⁺² content. Conversely, the y-axis separates the waters of the first group (Al-Hujrah zone) with negative values, indicating that these waters have a higher content of SO₄⁻² compared to Cl⁻. In contrast, the waters of groups (2) and (3) exhibit higher concentrations of Cl⁻⁻ than SO₄⁻² (downstream part of Al-Hujrah and Ranyah zone).

4.3 Natural processes

Various techniques, including the Gibbs and molar ratio diagrams with a common denominator are commonly used to study the diversity in water composition. However, these diagrams, which rely on molar ratios for their coordinates, cannot be effectively analyzed using statistical methods due to their lack of orthogonal directions. This limitation hinders the modeling stage. To better understand the natural processes controlling the chemistry of groundwater, the revised Gibbs diagram was employed in this study. This version of the Gibbs diagram is consistent with the kind of compositional data recently introduced by (Buccianti, 2015) through the use of isometric log-ratio equilibria for groundwater samples.

The groundwater geochemical analysis of Wadi Ranyah has been presented using a binary diagram. The coordinates in this diagram are determined by the equations listed in Supplementary Materials (Table S2). Fig. 5 depicts the three Gibbs geochemical mechanisms, as well as the position of the seawater composition. In this diagram, lower values on the y-axis indicate a low concentration of elements in the aqueous solution, while higher values indicate the prevalence of water–rock interaction processes. Additionally, the x-axis demonstrates the dominance of Na⁺ over Ca⁺², as illustrated in Fig. 5. The majority of water samples are in the water–rock interaction domain, with two distinct groups. The first group consists of wells located in the upstream part (Al-Hujrah zone). These wells are characterized by a dominance of Ca⁺² over Na⁺, indicating the influence of water-aquifer reservoir interaction. The second group includes wells from the downstream part (Ranyah zone), where most points lie along a vertical line representing an equilibrium between Na⁺ and Ca⁺². In this group, the y-axis values continue to increase, suggesting that water chemistry is affected by evaporation. Wells H51 and H55 are located close to the evaporation line, where an increase in Ca⁺² is generally observed due to anthropogenic activities.

Close

Fig. 5

Isometric log-ratio balances for Gibbs's diagram revised as proposed by Buccianti (2015).

Export to PPT

4.4 Water quality index

In the present study, we calculated the groundwater quality index for domestic use in the Wadi Ranyah area. We used weights based on relative loads and eigenvalues for trace metal elements and major ions, which were converted to central logarithmic ratio (clr) as input variables in the principal components (PC) analysis. The standards set by (WHO, 2017) have also been taken into account, as shown in Table 3. The water quality index (WQI) tool we employed shows great promise. It can be employed to evaluate overall water quality or serve as a pollution index. Furthermore, this tool facilitates swift and cost-effective monitoring of multiple water bodies, helping identify sources of pollution or elements that require immediate monitoring or treatment (Ewuzie et al., 2021).

In this research, the estimated water quality index (WQI) values (Supplementary Materials Table S6) for Group (1) ranged from 1410 to 1948, with an average of 1731, while Group (2) had a range of 1260 to 1930, with an average of 1502 (Fig. 6). Group (3) showed a WQI range of 1412 to 1957, with a mean of 1712. It was observed that all groundwater samples from Wadi Ranyah were categorized as water of very poor quality and unsuitable for consumption (WQI > 300). On the basis of the sub-indices, it was noted that lead had high sub-index values (>24). Consequently, the high WQI values correspond to higher levels of this parameter in groundwater, where lead concentrations range from 0.24 to 0.4 mg/L, well above the WHO limit (0.01 mg/L).

Close

Fig. 6

Graphic representation of the water quality index in the study area.

Export to PPT

A significant correlation was found between WQI and fluorine and manganese (0.45 > R² > 0.69), while a strong correlation was observed with lead (R² = 0.99) (Table 4). This indicates that water quality is influenced only by high concentrations of pollution elements (F, Mn and Pb). These elements can originate from two main sources: (i) weathering of crystalline and metamorphic rocks such as granites and gneisses (Subyani and Al Ahmadi, 2010), and (ii) human activities such as road traffic and agriculture. To guarantee human health, it is essential to regularly monitor groundwater quality, implement water treatment measures where necessary, and provide the populations concerned with alternative sources of clean, safe drinking water.