Introduction

As an archipelago surrounded by volcanoes, volcanic aquifers play an important role in meeting the water needs of Indonesia (Toulier et al. 2019; Maria et al. 2021; Fadillah et al. 2023). The potential for groundwater in volcanic aquifers in tropical countries is significant because it is the primary source of aquifer recharge due to high rainfall intensity (Razi et al. 2023). In addition to its abundant availability, groundwater quality is often better than surface water because, during the recharge to abstraction process, it is not significantly affected by anthropogenic activities on the surface (Sunkari et al. 2022). Groundwater resources in the Yogyakarta-Sleman Groundwater Basin are widely used to meet people’s water needs, particularly for domestic, industrial, and agricultural activities. The daily groundwater abstraction in this area reaches 340,000 m3/day (Hendrayana et al. 2023). However, the massive utilization of groundwater can lead to several environmental impacts, such as the reduction of groundwater reserves, land subsidence, drought in the dry season, and groundwater pollution (Barua et al. 2020; Wilopo et al. 2021; Taftazani et al. 2022).

The complexity of volcanic aquifer systems has encouraged comprehensive studies on aquifer characteristics, groundwater flow systems, and geochemical processes (Saberinasr et al. 2019; Selmane et al. 2022; Patel et al. 2023). Hydrogeochemical studies in volcanic aquifers are important for providing an overview of the interaction between groundwater and aquifer minerals and the groundwater circulation process from recharge to discharge zones (Dragon 2021). Previous hydrogeochemical studies have demonstrated that as groundwater travels from recharge to discharge zones, its chemical composition and overall quality can be modified by various geogenic factors (Dong and Gao 2022; Liu et al. 2022). Groundwater is naturally in contact with rocks and soils, allowing for the presence of various minerals in both dissolved ion and non-ionic solutions (Ansari et al. 2019). As a result, water–rock interaction processes such as mineral dissolution, weathering, evaporation, and ion exchange also shape groundwater geochemistry. (Liu et al. 2022; Sunkari et al. 2022; Ribinu et al. 2023). A deeper understanding of the water–rock interaction process can be obtained by inverse modeling to simulate the geochemical composition of groundwater and determine the main processes that affect groundwater quality in the study area (Rafiq et al. 2022). In addition, the geochemical content of groundwater also depends on anthropogenic activities at the surface, climatic conditions, and the quality of recharge water (Pant et al. 2018; Ansari et al. 2019; Wilopo et al. 2021).

Research related to hydrogeochemistry and groundwater quality has been conducted several times by previous researchers in the Yogyakarta-Sleman Groundwater Basin, and it has continued to develop over the past few decades. Early research was conducted by MacDonald and Partners (1984), who examined the general hydrogeology of the Yogyakarta-Sleman Groundwater Basin, followed by the analysis of aquifer interconnections in the Merapi Spring Belt (Boulom et al. 2015) and hydrogeochemical characteristics and groundwater quality on Merapi’s slope (Hendrayana et al. 2023). Most of these studies were conducted locally using limited analytical parameters. Therefore, this study comprehensively investigated the hydrogeochemical characteristics, groundwater mineralization processes, and hydrogeochemical evolution along groundwater flow paths on a basin scale. Specifically, this study aims to characterize the hydrogeochemical evolution of multilayer aquifer systems, reveal the dominant hydrogeochemical processes and mineral sources of groundwater, and understand aquifer interconnections in groundwater flow systems in unconfined and confined aquifers. The results of this study will help interpret the hydrogeochemical evolution of complex multilayer volcanic aquifer systems. Integrating this information is part of scientific-based evidence to achieve Sustainable Development Goals (SDGs) 6 (WWAP 2018).

Study area

Geographical setting

The research location follows the boundary of the Yogyakarta-Sleman Groundwater Basin, which includes the southern slope of Mt. Merapi (Fig. 1a). This study area includes three administrative zones, namely Sleman Regency, Yogyakarta City, and Bantul Regency, all of which are part of the Special Region of Yogyakarta Province, Indonesia. The Yogyakarta-Sleman Groundwater Basin is bounded by two major rivers: the Opak River in the east and the Progo River in the west, and the Indian Ocean directly borders the southern part (Hendrayana and Vicente 2013). The groundwater level in this area generally follows the season, where during the rainy season from January to April, the groundwater reaches the shallowest level and then begins to decline gradually when entering the dry season (Razi et al. 2023). Furthermore, the groundwater table spread radially and corresponded to the morphological pattern of the volcano. Towards the south, there is a decrease in the topographic gradient and a reduction in the hydraulic gradient. It causes the groundwater flow velocity to decrease towards the south of the Groundwater Basin (Putra and Indrawan 2012).

Fig. 1
figure 1

a Map showing groundwater sampling sites from wells, springs, and artesian wells in the Yogyakarta-Sleman Groundwater Basin, b North–south lithological cross-section from integrated monitoring wells

Full size image

Geology and hydrogeology of study area

The Yogyakarta-Sleman Groundwater Basin area generally comprises volcanic deposits from several main formations: the Old Merapi Deposit, Sleman Formation, and Yogyakarta Formation (MacDonald and Partners 1984). The Old Merapi Deposit consists of basalt and andesite lavas with indurated breccias exposed around the upper cone of Mt. Merapi and deposited during the Upper Pleistocene. Meanwhile, the Sleman Formation was estimated to be Upper Pleistocene to Holocene in age (MacDonald and Partners 1984). The Sleman Formation has been determined to be the lower part of a volcaniclastic unit previously included in the Young Merapi Volcanic Formation. To the north of the upper slopes of Merapi, it consists of sand and gravel, with boulders derived from volcanic outburst products. Towards the south, the Sleman Formation is overlain by the Yogyakarta Formation, which consists of interbedded sand, gravel, silt, and clay sequences. The basement of the basin is dominated by limestone, consisting of marl, tuff, and conglomerate as part of the Sentolo Formation.

The hydrogeological system formed by the Yogyakarta Formation and Sleman Formation in the Yogyakarta-Sleman Groundwater Basin is called the Merapi Aquifer System and includes a multilayered aquifer (Ilham et al. 2021). These layered aquifers have relatively similar hydraulic properties and are interconnected. The Yogyakarta-Sleman Groundwater Basin’s aquifers are mostly grouped into upper (unconfined) and lower (semi-confined) aquifers. The hydraulic conductivity in the Yogyakarta-Sleman Groundwater Basin ranges from 0.8 to 95 m/day, with an average value of about 8.6 m/day. However, the specific yield of the upper aquifer is estimated to be approximately 20%, which can be classified as a productive aquifer (Hendrayana et al. 2023). The hypothetical recharge area is located between elevations of 2968–700 m above sea level (masl), the transition area between elevations 700–200 masl, and the discharge area has an elevation between 200 to 0 masl (Hendrayana and Vicente 2013).

Methodology

Hydrogeochemical sampling and analysis

Groundwater samples were obtained from 35 observation stations representing wells, springs, and artesian wells. There were 23 wells representing the dug wells and shallow boreholes, eight springs, and four artesian wells (representing a confined aquifer). Field data were collected considering seasonal variability, with 66 rainy and dry season samples. Groundwater sampling is very susceptible to contamination; therefore, the sampling process must be performed using airtight, high-density linear polyethylene (HDPE) plastic bottles. Each sample taken at each location was placed into a cool box that served as a storage container to maintain the sample quality until it was analyzed in the laboratory. The samples were preserved with HNO3 until a pH of < 2 for cation analysis. Before further laboratory analysis, groundwater samples were filtered using a 0.20 μm filter membrane to separate dissolved substances and impurities in the samples. The physical parameters observed included temperature (T), acidity (pH), electrical conductivity (EC), and total dissolved solids (TDS), which were measured in situ with the HANNA device (HI 9813-6). Furthermore, groundwater samples were analyzed using 850 Metrohm Professional Ion Chromatography (IC) to obtain the major ions of the groundwater, such as Ca, Mg, Na, K, SO4, Cl, and NO3. Simultaneously, HCO3 was calculated using the volumetric titration method. The percentage charge balance error (% CBE) was calculated using Eq. 1 to ensure the reliability and accuracy of the measured chemical parameters.

$$\%CBE= \frac{{N}_{c}-{N}_{a}}{{N}_{c}+{N}_{a}} \times 100$$
(1)

where \({N}_{c}\) and \({N}_{a}\) are the cation and anion concentrations (in meq/L).

Hydrogeochemical data were analyzed using Trilinear Piper and Chada diagrams to determine groundwater facies and the main chemical processes occurring in the Yogyakarta-Sleman Groundwater Basin. Most of the graphs used in this study were generated using Microsoft Excel and AquaChem version 4.0 software, and statistical processing was performed using IBM SPSS Statistics 25.0. Correlation analysis between parameters was conducted based on the Pearson correlation, which allows the grouping of groundwater based on the similarity of chemical and physical properties in groundwater samples (He et al. 2023). To determine the process of water–rock interaction, ion exchange, and its implications on the chemical content of groundwater, an analysis with a Gibbs diagram was conducted, strengthened by the interpretation of the bivariate plot of major groundwater ions. This approach has been widely used to understand groundwater’s hydrogeochemical evolution at the basin scale (Kouadra and Demdoum 2020).

Hydrogeochemical modelling

The saturation index (SI) was calculated using PHREEQC 3.0, a geochemical modeling tool for evaluating groundwater thermodynamics and describing the saturation index equilibrium of any given mineral (Ansari et al. 2019). The SI of various materials is controlled by the ion activity product (IAP) and solubility product (K). The relationship between the SI and different minerals is defined as the ratio of IAP to K:

$${\text{SI}}={\text{log}}\left(\frac{IAP}{K}\right)$$
(2)

The SI results express the chemical equilibrium trends of water and minerals and the interaction process between groundwater and rock. If unsaturated (SI < 0), minerals will continuously weather in the groundwater; if oversaturated (SI > 0), minerals will precipitate; and if SI is close to 0, minerals are in a fixed phase in an equilibrium state (Gao et al. 2019).

Result and discussion

Physicochemical characteristics

The physicochemical characteristics of the groundwater for each hydrogeochemical parameter representing the unconfined and confined aquifers are shown in Table 1. TDS is an important parameter that represents the dissolved minerals in groundwater and can provide an initial condition for the mineralization process. The TDS values for the unconfined aquifer ranged from 50 to 320 mg/L, with a mean value of 159.84 mg/L. In comparison, the TDS value for the confined aquifer ranged from 320 to 480 mg/L (mean 395 mg/L). The higher TDS value for groundwater in the confined aquifer (represented by artesian wells) indicates a more intensive water–rock interaction process compared to groundwater in an unconfined aquifer. There was also a tendency for TDS values to increase along the flow path from the recharge to the discharge area. The same pattern also occurred for EC. For the unconfined aquifer, the EC value ranged from 120 to 660 µS/cm (mean 335.41 µS/cm), whereas for the confined aquifer, the value ranged from 320 to 480 µS/cm (mean 395 µS/cm). Low TDS and EC values were obtained from shallow dug well samples that received direct rainwater recharge. This process is controlled by the direct infiltration of meteoric water into the unsaturated layer; therefore, the contact time between the groundwater and rock is relatively short (Khettouch et al. 2022). The pH value for the unconfined aquifer showed variations between 5.7 and 7.70 (mean 6.59), whereas for the confined aquifer, the value was greater and ranged from 7 to 8. The high pH value for deep artesian wells (total depth 90 to 130 m) indicates contact between groundwater and the basement of the basin, which is composed of Tertiary limestone.

Table 1 The statistics of physicochemical measurements on groundwater samples
Full size table

The geochemical content of groundwater under natural conditions is controlled by the main cations of Ca, Mg, Na, and K and the main anions of HCO3, SO4, and Cl, whereas NO3 is more commonly derived from anthropogenic activities on the surface (such as agriculture and domestic waste) (Gao et al. 2019; Saha et al. 2020; Singh et al. 2020). In general, the average groundwater cation concentrations for unconfined aquifers followed the order Ca > Na > K > Mg, while for confined aquifers followed the order Na > Ca > Mg > K. There was no significant difference in cation concentration for the rainy and dry season both for unconfined and confined aquifers as shown in Fig. 2. Ca and Na were the main cations found in groundwater samples. The concentration of Na in an unconfined aquifer ranged from 3.76 to 64.82 mg/L (mean 21.34). In contrast, the Na concentration for confined aquifers ranged from 121.75 to 155.57 mg/L. The high Na content in deep artesian wells indicates silicate weathering and water’s interaction with aquifer materials (Saha et al. 2019). Ca concentration for unconfined aquifers ranged from 10.99 to 62.14 mg/L (mean 32.87 mg/L), and for confined aquifers, a higher Ca value ranged from 20.61 to 65.59 (mean 37.26 mg/L). The concentration of Mg for unconfined aquifers ranged from 2.03 to 20.29 mg/L (mean 9.63 mg/L), while for confined aquifers ranged from 8.47 to 27.68 mg/L (mean 16.20 mg/L). Reactions involving Mg dissolution are influenced by the amount of CO2 in groundwater under certain conditions derived from the decomposition of ferromagnesian minerals, such as olivine, pyroxene, and amphibole (Saha et al. 2019). K is one of the natural elements in the main primary ions of groundwater, but its concentration is relatively lower than that of Ca, Mg, and Na. The K concentration for the unconfined aquifer varied from 1.86 to 47.07 mg/L (mean 11.89 mg/L), and in the confined aquifer varied from 7.46 to 20.85 mg/L (mean 13.77 mg/L).

Fig. 2
figure 2

Boxplot diagram showing the variability of major ions in groundwater during the dry and rainy seasons

Full size image

The average anion concentration in the unconfined aquifer followed the order HCO3 > SO4 > NO3 > Cl. Meanwhile, the anion concentration for confined aquifer followed the order HCO3 > Cl > SO4 > NO3. HCO3 was the main anion in both unconfined and confined aquifers. As shown in Fig. 2, there was no significant change in anion concentration between the rainy and dry seasons, except for HCO3. Various geochemical processes, such as carbonate dissolution from the atmosphere and weathering of silicate minerals, are the dominant factors causing high HCO3 concentrations in groundwater (Ribinu et al. 2023). The enrichment of HCO3 over Cl during the rainy season reveals that groundwater chemistry is influenced by the recharge process and weathering of carbonate and silicate minerals (Manikandan et al. 2020). A previous study also reported that flushing weathered materials from the vadose zone during the rainy season resulted in a higher concentration of HCO3 in groundwater than other major ions (Rao 2006). HCO3 concentration for the unconfined aquifer ranged from 48 to 375 mg/L (mean 148.20 mg/L); meanwhile, for the confined aquifer, HCO3 ranged from 353.80 to 439.20 mg/L (mean 387.35 mg/L). Cl concentration for unconfined aquifer varied from 1.77 to 62.08 mg/L, with a mean value of 14.48 mg/L. Meanwhile, the Cl value was relatively higher in the confined aquifer, ranging from 17.58 to 147.79 mg/L (mean 96.75 mg/L). Hence, it can be proven that the Cl concentration tends to increase with regional groundwater circulation. Another important major ion in the groundwater is SO4. High SO4 concentrations may be due to dissolved halite, anhydrite, and gypsum minerals in groundwater (Ghani et al. 2022). For an unconfined aquifer, SO4 concentration varied from 7.85 to 85.90 mg/L (mean 26.97 mg/L), while the SO4 concentration from 3.22 to 107.48 mg/L (mean 41.51 mg/L). However, NO3 was the lowest anion concentration in a confined aquifer, ranging from 0.05 to 2.27 mg/L (mean 0.70 mg/L), and the second lowest anion in an unconfined aquifer ranged from 0 to 85.96 mg/L (mean 16.73 mg/L).

The interrelationships between the 11 physiochemical parameters were evaluated using Pearson’s correlation. This correlation matrix is widely used to estimate the degree of dependency of a variable on another. As shown in Table 2, TDS and EC showed highly significant correlations with several groundwater chemical parameters, such as Na, Cl, HCO3, Mg, and Ca. There was also a significant correlation between Na and Cl (0.802), Na and HCO3 (0.820), Ca and Mg (0.788), and Mg and HCO3 (0.746), indicating that many ions in the groundwater may originate from halite, gypsum, and carbonate dissolution (He et al. 2023; Nsabimana and Li 2023). All major ions were low and negatively correlated with NO3. High NO3 concentrations were identified in shallow dug wells in the Yogyakarta City agglomeration area. Previous studies have revealed that the nitrate concentration in the Yogyakarta-Sleman Groundwater Basin mainly comes from domestic waste due to improper sewage systems (Fathmawati et al. 2017; Putra 2015). It suggests that NO3 in groundwater is not derived from natural mineralization processes but is influenced by anthropogenic activities at the surface, especially in areas with high population density.

Table 2 Pearson correlation analysis for each physiochemical parameter in dry and rainy seasons (n = 66)
Full size table

Groundwater facies

A Trilinear Piper diagram characterized the groundwater facies and geochemical processes. The Trilinear Piper diagram is widely used to determine water type or hydrogeochemical facies by plotting major ions on subdivisions of diamond-shaped planes. Therefore, this concept is widely used to differentiate naturally occurring water and rock interactions (Davraz and Batur 2021). Based on the Trilinear Piper diagram (Fig. 3a), groundwater facies in the Yogyakarta-Sleman Groundwater Basin are Ca–Mg–HCO3, mixed Ca–Na–HCO3, mixed Ca–Na–Cl, Na–HCO3, and Na–Cl–SO4 type. About 75% of the groundwater in the study area was confirmed to be Ca–Mg–HCO3 type. This type is dominated by shallow groundwater in unconfined aquifers with high HCO3 and relatively low TDS. This condition is consistent with previous research in the southern part of Merapi, which revealed that shallow groundwater on the slopes of Merapi is dominated by HCO3, Ca, and Mg concentrations (Boulom et al. 2015; Hendrayana et al. 2023). HCO3 is present in limestone, dolomite, and secondary forms, filling cracks in aluminosilicate rocks in the form of cement and crusts in soils (Clark 2015). The meteoric water permeating the aquifer layer is enriched by CO2, which may be retained by the dissolution of CO2 from the atmosphere or the degradation of plant debris (Delkhahi et al. 2020). When CO2 dissolves in groundwater, it forms carbonic acid, which naturally serves as the major acid for dissolving carbonates or aluminosilicate rocks (Mitchell et al. 2010). In this case, HCO3, Ca, and Mg are released into groundwater during the infiltration process of meteoric water (Eqs. 3 and 4):

Fig. 3
figure 3

Groundwater facies classifications using a Trilinear Piper diagram and b Chadha diagram for groundwater facies classification

Full size image
$${{\text{CaCO}}}_{3}+ {{\text{H}}}_{2}{{\text{CO}}}_{3} \to {{\text{Ca}}}^{2+} + 2{{\text{HCO}}}_{3}^{-}$$
(3)
$${{\text{CaMg}}({{\text{CO}}}_{3})}_{2}+2{{\text{H}}}_{2}{{\text{CO}}}_{3}\to {{\text{Ca}}}^{2+}+{{\text{Mg}}}^{2+}+4{{\text{HCO}}}_{3}^{-}$$
(4)

The extent of mineralization is a function of the \({P}_{{CO}_{2}}\) of soils in the recharge area and the degree of openness of the system during weathering (Clark 2015; Delkhahi et al. 2020).

Groundwater in shallow aquifers evolved into a mixed Ca–Mg–Cl type towards discharge zones, representing meteoric water that has evolved during the relatively long circulation (Nsabimana and Li 2023). In the confined aquifer, there was a significant increase in Cl, Na, and SO4 ions due to the influence of the water–rock interaction. The groundwater facies changed from mixed Ca–Na–HCO3 to Na–HCO3 and Na–Cl–SO4. The enrichment of Na content is presumably due to interaction with fluvio-volcanic deposits from the Merapi Aquifer System rich in plagioclase and hornblende, as well as the ion exchange process (Boulom et al. 2015). The natural increase in Cl as a conservative ion indicates the salinization process for regional groundwater flow in the lower aquifer (Huang and Ma 2019). This process is caused by the weathering of aquifer rock minerals in some artesian wells (SA-1, SA-2, SA-3, and SA-4).

In addition to using Trilinear Piper diagrams, this study used a classification based on the Chada diagrams introduced by Chadha (1999). The Chadha diagram is an updated Trilinear Piper diagram and an extension of the Durov diagram. In addition to providing a classification of groundwater facies, the Chada diagram can show natural groundwater geochemical processes, including the directional exchange process, recharge water, saline water, and inverse ion exchange processes. As shown in Fig. 3b, most groundwater samples in the rainy and dry seasons had the Ca–Mg–HCO3 type derived from recharge water. All groundwater in these facies is shallow groundwater that belongs to the unconfined aquifer in the Yogyakarta-Sleman Groundwater Basin system. However, there is also some shallow groundwater with Ca–Mg–Cl facies in the discharge area, where Cl replaced the dominance of HCO3. Meanwhile, groundwater originating from confined aquifers is characterized by Na–HCO3 facies, and one artesian well has Na–Cl–SO4 facies. The ion exchange process controls groundwater in the lower aquifer.

Rock weathering and mineral dissolution

The Gibbs diagram can be used to infer the main geochemical processes (Pant et al. 2018). This method can reveal geochemical evolution involving various processes, such as precipitation, rock weathering, and evaporation-crystallization processes. The Gibbs diagram is widely used to give the initial condition of the main geochemistry processes for further analysis and interpretation (Verma and Singh 2021; Mohamed et al. 2022). The Gibbs diagram is constructed based on the equivalence of Na/(Na + Ca) and Cl/(Cl + HCO3) ratios versus TDS (Zhang et al. 2020).

$$\text{Gibbs ratio I }(\text{for cation})=\left[\frac{{{\text{Na}}}^{+}}{{({\text{Na}}}^{+}+ {{\text{Ca}}}^{2+})}\right]$$
(5)
$$\text{Gibbs ratio II }(\text{for anion})=\left[\frac{{{\text{Cl}}}^{-}}{({{\text{Cl}}}^{-}+ {{\text{HCO}}}_{3}^{-})}\right]$$
(6)

In the Gibbs diagram, if a groundwater sample is plotted at the bottom right with a low TDS value but higher Na/(Na + Ca) and Cl/(Cl + HCO3) values, it indicates that the chemical composition is mainly controlled by atmospheric precipitation and undergoes a dominant weathering process. In contrast, a high TDS and high Na/(Na + Ca) and Cl/(Cl + HCO3) values indicate evaporation’s dominant influence. As the Gibbs diagram (Fig. 4a) reveals, water–rock interaction is the primary process affecting the geochemical composition of groundwater in the Yogyakarta-Sleman Groundwater Basin for samples from spring, well, and artesian wells. This trend occurs in both the rainy and dry seasons. The rock type strongly influences the chemical concentration of groundwater in the aquifer during the residence time (Ansari et al. 2019). Further interpretation was performed using Gaillardet’s diagram to confirm the water–rock interaction process in groundwater. The diagram classifies groundwater interaction processes as controlled by evaporite dissolution, silicate weathering, and carbonate dissolution (Gaillardet et al. 1999). Based on Fig. 4b, all samples fall into the silicate weathering section, indicating that the dissolution of silicate minerals plays a significant role in groundwater geochemistry in the Yogyakarta-Sleman Groundwater Basin. Albite is a mineral that generally has a weathering process in groundwater with the following reaction (Zhang et al. 2020):

$${\text{2Na}}^{ + } \left[ {{\text{AlSi}}_{{3}} {\text{O}}_{{8}} + {\text{ 2CO}}_{{2}} } \right]^{ - } + {\text{11H}}_{{2}} {\text{O}} \to {\text{ 2Na}}^{ + } + {\text{2HCO}}_{{3}}^{ - } + {\text{Al}}_{{2}} {\text{Si}}_{{2}} {\text{O}}_{{5}} \left( {{\text{OH}}} \right)_{{4}}^{{}} + {\text{4H}}_{{4}} {\text{SiO}}_{{4}}$$
(7)
Fig. 4
figure 4

a Gibbs diagram showing the mechanism controlling groundwater chemistry in the study area from spring, well, and artesian well samples. b Gaillardet diagram using Na-normalized molar ratios in the dissolved phase

Full size image

The abundance of major ions is influenced by silicate weathering or large amounts of CO2 dissolved in meteoric water. In this case, the weathering of carbonates produces Ca, Mg, and HCO3, while the weathering of silicate minerals produces Ca, Mg, Na, K, Si, and HCO3. At the same time, the dissolution of evaporites can produce Ca, Mg, Na, K, SO4, and Cl (Pant et al. 2018). However, Gaillardet’s diagram requires a more comprehensive classification technique of water–rock interaction processes, as shown in Fig. 5.

Fig. 5
figure 5

Scatter plot showing the correlation between major ions and hydrogeochemical processes. a Ca vs. SO4, b Na + K vs. Cl + SO4, c Cl vs. Na, and d total cation vs. (Na + K)

Full size image

Various approaches can be used to compare the ratios of major ions using bivariate plots, such as the comparison of Ca vs. SO4, Na + K vs. Cl + SO4, Cl vs. Na, and total cations vs. (Na + K). A comparison between Ca and SO4 (Fig. 5a) reveals that almost all groundwater samples in the study area fell below the 1:1 line, indicating that silicate weathering was the primary process controlling the chemical content of the groundwater. If Ca and SO4 originate from gypsum, the ratio is 1:1 (Zhang et al. 2020). However, only three samples were found to originate from the gypsum/anhydrite or sulfate dissolution process. The ratio of Cl + SO4 to Na + K was calculated to confirm the silicate weathering process (Fig. 5b). If the groundwater samples fall below the 1:1 line, it can be proved that the geochemistry of the groundwater is highly dependent on the silicate weathering process.

Figure 5.c shows the molar ratio between Cl and Na, where a Cl/Na ratio close to 1 indicates that halite dissolution contributes to the Na concentration in the groundwater. Meanwhile, if the ratio exceeds 1:1, silicate weathering dominates the geochemical process (Ansari et al. 2019). The contribution of cations to silicate weathering in groundwater can also be estimated using the ratio of (Na + K) to total cations (Fig. 5d). In the research area, groundwater samples in the unconfined aquifer fall below the 1:1 line, while artesian well samples are above the 1:1 line. The low (Na + K) concentration may be due to the Ca/Na exchange, which can reduce the Na concentration in the groundwater (Lakshmanan et al. 2003). Four scatter plots in Fig. 5a–d support the interpretation that silicate weathering is the primary geochemical process affecting groundwater. Simultaneously, carbonate dissolution did not significantly contribute to the presence of major ions.

Ion exchange process

Ion exchange is one of the main processes in groundwater systems that influence hydrogeochemical evolution (Saha et al. 2020). Sedimentary mineral surfaces in aquifers, particularly the clay fraction, carry a negative charge that can absorb water-soluble cations. The chemical process changes the equilibrium between the adsorbent and solution phases, which causes cations to be adsorbed into or replaced by new cations in the groundwater (Carol et al. 2012). The alteration of the chemical content of groundwater can cause a shift in the equilibrium between the solution and solid phases, leading to the adsorption of cations into new positions or the replacement of cations in the groundwater. Based on Eq. 8, the ion exchange process can occur when calcium-rich water flows through a sodium-loaded clay aquifer, where Ca or Mg replaces Na during the interaction between groundwater and minerals. This process transforms groundwater from a Ca–Mg–HCO3 type in the recharge area to a Na–HCO3 type in the discharge area, which causes the position of the samples on the Piper diagram (Fig. 3.a) to shift gradually to the lower right and middle part. Meanwhile, the reverse ion exchange process follows Eq. 9, where Na exchanges Ca (or Mg) minerals with minerals in the groundwater aquifer.

Ion exchange:

$${{\text{Ca}}}^{2+}\left({{\text{Mg}}}^{2+}\right)+2{\text{Na}}-{\text{X}}\to 2{{\text{Na}}}^{+}+{\text{Ca}}\left({\text{Mg}}\right)-{{\text{X}}}_{2}$$
(8)

Reverse ion exchange:

$$2{{\text{Na}}}^{+}+{\text{Ca}}\left({\text{Mg}}\right)-{{\text{X}}}_{2} \to {{\text{Ca}}}^{2+}\left({{\text{Mg}}}^{2+}\right)+2{\text{Na}}-{\text{X}}$$
(9)

The cation exchange process was determined by plotting the ratio of (Ca + Mg)/(Na + K) vs. total cations (in meq/L). As shown in Fig. 6a, most of the groundwater in the unconfined aquifer experienced a reverse ion exchange process. In contrast, groundwater samples in the confined aquifer are confirmed to have experienced the ion exchange process. Furthermore, Fig. 6b shows the bivariate plot between Na vs. Ca + Mg (in meq/L), revealing that the alkaline earth element dominates over alkali (Zaidi et al. 2015). If the sample is on the right side of the line (1:1), it indicates forward ion exchange, while on the left, it indicates reverse ion exchange. From Fig. 6b, only four samples show an ion exchange trend, dominated by groundwater in confined aquifers that contain relatively high concentrations of Na + K compared to other ions.

Fig. 6
figure 6

Relationship among chemical parameters in groundwater that explain ion exchange and reverse ion exchange process. Relation of the a sum of cation vs. (Ca + Mg)/(Na + K), b Na vs. (Ca + Mg), c (Ca + Mg) vs. (HCO3 + SO4), d (Na – Cl) vs. (Ca + Mg – (HCO3 + SO4))

Full size image

The ionic interaction process in the aquifer can also be recognized using the ratio of Ca + Mg vs. HCO3 + SO4, as shown in Fig. 6c. Groundwater samples along the 1:1 line represent the dissolution of calcite, gypsum, and dolomite, while samples falling below the 1:1 line indicate the reverse ion exchange process. During the ion exchange process, Ca is retained in the aquifer while Na is released to groundwater, whereas in the reverse ion exchange process, Na is retained in the aquifer while Ca is released to groundwater (Mahmoudi et al. 2017). The process is confirmed by Fig. 6d, which shows the bivariate plot between Na − Cl vs. Ca + Mg − (HCO3 + SO4). All groundwater samples in the study area tend to release Na and adsorption of Ca. The upper left sector has a deficit of Na compared to Ca and Mg, as shown in Fig. 6d. This sector’s groundwater is highly mineralized, but sodium is lower than chloride. The lower right sector shows an excess of Na concerning Cl and a deficiency of Ca and Mg compared to the upper right section (Gning et al. 2017). Therefore, understanding the ion exchange and reverse ion exchange processes can provide comprehensive information about the hydrogeochemical evolution of groundwater from the recharge to the discharge zones.

Saturation index

Hydrogeochemical processes along the groundwater flow direction are primarily determined by mineral weathering over long residence times (Kouadra and Demdoum 2020). Calculating the mineral equilibrium can indicate the thermodynamic processes occurring in groundwater, which can be determined by the SI values (Gao et al. 2019). The SI value was obtained based on Eq. 2, where SI < 0 indicates undersaturated, SI = 0 indicates equilibrium, and SI > 0 indicates oversaturated. An undersaturated state suggests that the mineral has the potential to continue releasing ions into groundwater (Ghani et al. 2022). SI calculations were performed on several minerals: anhydrite, aragonite, calcite, dolomite, halite, and sylvite. In general, only aragonite (CaCO3), calcite (CaCO3), and dolomite (CaMg(CO3)2) were saturated in the groundwater of specific samples. The dissolution processes of these minerals can be expressed in the following reactions:

$${{\text{CaCO}}}_{3}\to {{\text{Ca}}}^{2+}+{{\text{CO}}}_{3}^{2-}$$
(10)
$${{\text{CaMg}}{({\text{CO}}}_{3})}_{2}\to {{\text{Ca}}}^{2+}+{{\text{Mg}}}^{2+}+ {2{\text{CO}}}_{3}^{2-}$$
(11)

Table 3 shows the maximum, minimum, mean, and standard deviation values for each groundwater sample. Saturated minerals only occurred in confined aquifers for aragonite, calcite, and dolomite. The dissolution of ions from aragonite, calcite, and dolomite minerals into groundwater in the Merapi Aquifer System was influenced by the relatively long groundwater flow and residence time, as well as the possibility of contact with carbonate basement rocks, specifically limestone from the Sentolo Formation. In contrast, no samples from the unconfined aquifer were saturated with the following minerals during either the rainy or dry seasons.

Table 3 Statistical summary of SI of certain minerals for groundwater samples
Full size table

Figure 7a shows the relationship between SI and TDS for the most prevalent carbonate minerals: dolomite, calcite, gypsum, and anhydrite. It is commonly observed that TDS exhibits a positive correlation with SI values, as there is a tendency for SI values to escalate with an increase in TDS. Samples with TDS values < 300 mg/L were found to be mostly unsaturated with any minerals; therefore, minerals continuously weather in the groundwater at these locations. Figure 7b presents the relationship between the SI of dolomite vs. the SI of calcite, the common carbonate minerals found in groundwater. Most groundwater fell in dolomite and calcite unsaturated sections, whereas only artesian wells were identified in calcite and dolomite saturated conditions. The decrease in the concentrations of Ca and Mg due to the precipitation of dolomite and calcite resulted in a significant amount of free CO2 (Wu et al. 2021). The excessive dissolution of CO2 in groundwater leads to the dissolution of Na and K feldspars, which can be found in some artesian wells (Manikandan et al. 2020). This finding helps to elucidate the conditions of hydrogeochemical evolution in multilayer aquifers, where the groundwater originating from direct recharge has a lower mineral content, which then infiltrates and follows the regional flow system until it has sufficient time to reach a state of equilibrium and saturation in the aquifer.

Fig. 7
figure 7

SI graph a TDS vs. SI graph for selected minerals (dolomite, calcite, aragonite, and gypsum), b SI of calcite vs SI of dolomite

Full size image

Conceptual model

The hydrogeochemical characteristics of the study area explain the evolution of groundwater along the flow path from the recharge to discharge zones. In general, the geochemical evolution of groundwater from the recharge area to the discharge area gradually changes (Mahmoudi et al. 2017). Groundwater plays an active role in the current water cycle through short-distance flows to the discharge zone, where groundwater undergoes an active flushing process from rainwater through the porous rocks. Groundwater in the local flow system is characterized by dominant HCO3 with low TDS, as represented by the samples in the recharge zone (Fig. 8). Concurrently, the groundwater in regional flow systems must traverse long-distance aquifers composed of multilayer aquifers. Groundwater within the regional flow system is in a state of constant interaction with aquifer materials, which affects its chemical properties. In the transition and discharge zones, groundwater moves relatively slowly and produces dissolved solids, where SO4 and Cl become more dominant. This zone experienced a lower water flushing circulation than the recharge zone. The pressure in this zone is close to hydrostatic pressure, and the temperature tends to be constant (Dong and Gao 2022). Therefore, the groundwater gradually evolves from HCO3 type to SO4 type and then Cl type. Changes in anion species, as described previously, should be based on the scale and determination of physical and geological conditions.

Fig. 8
figure 8

Simplified conceptual model of hydrogeochemical evolution along the groundwater flow direction (a grey stiff diagram presents the sample from a confined aquifer, and a dark blue stiff diagram shows the sample from an unconfined aquifer)

Full size image

The existence of a hydraulic window from non-continuous impermeable layers in the aquifer system plays an important role in hydrogeochemical processes. Silt and clay are increasingly developed in the central area of the basin and can act as confining units between unconfined and confined aquifer systems (MacDonald and Partners 1984). Different hydraulic pressures between the lower aquifer (confined) and the upper aquifer (unconfined) cause water to flow from the lower aquifer to the upper aquifer in some parts of the basin (Putra 2015). This reflects the interconnection between the unconfined and confined aquifers in the Merapi Aquifer System. For confined aquifer, the high possibility of contact between groundwater and the basement rocks in this area, which is composed of limestone from the Sentolo Formation, can increase the ion dissolution from carbonate minerals, proven by the saturated conditions (SI > 1) for calcite, dolomite, and aragonite. These processes generate natural salinization, producing relatively high TDS and EC values.

Based on the hydrogeochemical evolution model (Fig. 8), the groundwater samples from the unconfined aquifer (MA-1, SB-2, SG-4, and SB-7) show a similar pattern of the stiff diagram, which indicates the similarity of the groundwater facies (Ca–Mg–HCO3). The hydrogeochemical evolution can be easily identified from four artesian wells representing the confined aquifer (SA-1, SA-2, SA-3, and SA-4). Artesian wells in the transition zone have mixed Ca–Na–HCO3 facies, indicating no significant change in groundwater facies but significant increases in pH, TDS, and EC as a consequence of the regional flow process from the recharge zone. The groundwater evolution in the confined aquifer continued in wells SA-2, SA-3, and SA-4, where the groundwater facies changed from mixed Ca–Na–HCO3 to Na–HCO3 and Na–Cl–SO4. The Mg in the confined aquifer replaced the absorbed Na in the aquifer medium. This phenomenon implies the occurrence of an ion exchange process. The rise in Cl ions in SA-4 at the end of the flow direction suggests a gradual alteration in the facies of water anions from HCO3 to SO4 and then to Cl, as demonstrated in Eq. 12. The geochemical evolution clarifies the groundwater circulation and the interconnection of the upper and lower aquifers in the Yogyakarta-Sleman Groundwater Basin.

Conclusion

In complex volcanic aquifer systems such as the Yogyakarta-Sleman Groundwater Basin, a comprehensive study of hydrogeochemical characteristics and groundwater circulation process is essential for sustainable groundwater management. This study provides new insights into multilayer volcanic aquifers’ conceptual geochemical evolution processes. The water–rock interaction processes in the regional groundwater flow can naturally cause groundwater salinization in the study area, resulting in relatively higher TDS and EC values. The increase in Cl ion at the end of the flow direction indicated a gradual change in the anion facies of the water from HCO3 to SO4 and then to Cl. Groundwater from the upper aquifer has dominant Ca-Mg-HCO3 facies, which are influenced by the direct infiltration process from meteoric water to the shallow aquifer. The groundwater facies evolved into a mixed Ca–Mg–Cl type along the flow direction towards the discharge zone. Simultaneously, when it reached the lower aquifer, the groundwater evolved into Ca–Na–HCO3, Na–HCO3, and Na–Cl–SO4 facies towards the discharge zone.

The Gibbs diagram demonstrates that water–rock interaction processes predominantly govern the variability of groundwater facies in the Yogyakarta-Sleman Groundwater Basin. It was confirmed by Gaillardet’s diagram, which demonstrated that all the samples were located in the silicate weathering section, indicating that silicate minerals are a major factor in groundwater geochemistry. Furthermore, the reverse ion exchange process dominates shallow groundwater in the recharge and transition zones. In contrast, the ion exchange process dominated the confined aquifer in the transition to the discharge zone. The Ca in groundwater is exchanged by Na, which is adsorbed in the clayey sediments associated with fluvio-volcaniclastic sediments from Merapi deposits. The mineralization of groundwater was also verified by the SI values calculated for certain minerals, with groundwater in the unconfined aquifer layer displaying an undersaturation state and groundwater samples in the confined aquifer showing saturation for aragonite, calcite, and dolomite minerals. The conceptual model generated from this study clearly explained the groundwater circulation throughout the basin and its physicochemical characteristics, which can provide scientific-based evidence for sustainable groundwater conservation. However, to confirm the origin and residence time during groundwater flow in the aquifer, future studies should combine groundwater hydrogeochemistry with tracer isotopes to obtain more comprehensive results.