The Control of Groundwater Flow Systems and Geochemical Processes on Groundwater Chemistry: A Case Study in Wushenzhao Basin, NW China

Abstract

The lowest reaches of a large-scale basin could be the discharge areas of local, intermediate and regional groundwater flow systems with significantly distinct travel distances and travel times. This study aims to delineate the groundwater chemical characteristics and the mechanism controlling the chemical evolution in the lowest reaches of the Wushenzhao Cretaceous basin, NW China. A total of 38 groundwater samples were collected and were chemically classified into five distinct water types by means of a Piper Plot. According to the hydrogeological setting and groundwater age, the spatial distribution of these water types is found to be associated with hierarchically nested groundwater flow systems (local and regional system): Types 1, 2, 3 and 4 belong to the local groundwater flow system, while type 5 belongs to the regional flow system. Graphical plots, stable isotopes and geochemical modeling techniques were used to interpret the observed compositions. The results show the dominance of carbonate and gypsum dissolution in type 1 waters; ion exchange in types 2, 3 and 4; and evaporite dissolution in type 5. In addition, human activities in the form of extensive irrigation also affect the chemical compositions of type 1 water. These findings are important for the sustainable management of groundwater resources in the study area.

1. Introduction

Geochemical processes taking place within the aquifer have a profound effect on groundwater chemical evolution [1,2]. An increased knowledge of groundwater flow and geochemical processes could contribute to the improved understanding of groundwater system in a region, which is of critical importance for the effective management and sustainable development of groundwater resources [3,4,5,6,7,8,9]. During the last decades, the integrated application of hydrochemistry and environmental isotopes has been proven to be a powerful tool for solving these hydrological problems [10,11,12,13].

Ordos Plateau is an important energy base located in the semiarid and arid regions, NW China (Figure 1a), which conducts plenty of industrial activities using substantial water resources [14]. Groundwater is the only reliable source of water for drinking, irrigation, industry and ecosystem [15]. In recent years, due to the rapid development of the energy base, the demand for water resources has increased tremendously. As a result, the groundwater level is increasingly declining, and the ecosystem environment is worsening, especially the Wushenzhao basin in the northeast of Plateau (Figure 1b) [16]. This may lead to the deterioration of groundwater quality in the future [8]. To help the sustainable development of groundwater (i.e., maintaining both yield and quality), a thorough understanding of the groundwater system, particularly the hydrochemical system, is essential in those Lake Basins.

The available geochemical investigations of groundwater have been just carried out at a regional scale in the Ordos Plateau [15,17,18], while detailed investigations on the groundwater geochemistry are rather scattered so far, and analyses conducted are limited to the major ions chemistry of groundwater in the Lake Basins, NE Ordos Plateau, such as Dake Lake, Mukai Lake and Hutongchahan (HTCH) lake (Figure 1b) [14,16,19,20,21,22]. In addition, the HTCH Lake Basin is chosen as the key study area (Figure 1c) because it is located in the lowest reaches of the Wushenzhao Basin and might serve as the fingerprint of nested flow systems. The purpose of this study is to delineate the flow and hydrochemical evolution of groundwater in the HTCH Lake Basin. A multi-tracer approach including ion chemistry and environmental isotopes (e.g., δ²H, δ¹⁸O, ⁸⁷Sr/⁸⁶Sr, δ³⁴S, ³H and ¹⁴C), together with geochemical modeling techniques, was employed. This work is the first detailed study on the groundwater chemistry in the study area and correlates the hydrochemical pattern with the groundwater flow systems. Moreover, groundwater age as well as improved potentiometric head contour is also first reported. These findings will significantly improve the sustainable development of groundwater in the study area.

2. Study Area

The Ordos Plateau is located in the north part of Ordos Basin, NW China (Figure 1a). It is one of the most important energy and industry bases due to an abundance of minerals resources, such as coal, natural gas, petroleum and halite [18,23]. The Plateau can be divided into three groundwater systems by the Sishi Ridge and Xinzhao Divide, namely the Wuding system, the Dosit system and the Molin system, as shown in Figure 1b [15]. The HTCH Lake Basin is located in the northeast part of Ordos Plateau and the north part of the Wuding groundwater system (Figure 1b). It belongs to the continental semiarid to arid climate: The annual precipitation averages about 340 mm/year; the mean annual potential evaporation is around 2300 mm/year, which far exceeds the precipitation.

Quaternary lacustrine sediments, aeolian sands and Cretaceous formations can be observed in the study area (Figure 2a). The first two are relatively thin, while the latter has a huge thickness of about 700–1000 m, which is further divided into the Huanhe and Luohe groups from the top to bottom (Figure 2a). The major aquifer is the poorly consolidated Cretaceous sandstone with a high permeability in the Plateau. Lens of mudstone are present in the Cretaceous formation, but they cannot constitute a regional aquitard due to their limited extent and thin thickness. Therefore, the Cretaceous aquifer is usually considered as a relatively homogenous and unconfined system [18]. The dominant minerals of aquifer rocks are quartz, Na- and K-feldspars, and calcite and dolomite [14].

The groundwater flow is driven by the wavy topography in the Ordos Plateau, which develops hierarchically nested groundwater flow systems. The local flow system generally develops between groundwater divides and neighboring depressions with a shallow circulation depth of around 200 m; intermediate flow systems occurs between the major divides and the relatively larger lakes with a relatively deep circulation depth; and the regional flow system has a much larger lateral extent and a higher circulation depth (up to 500 m) when compared with the local and intermediate flow systems [14,15,18,20,22].

There is a hypothesis of groundwater flow systems in the lake-concentrated regions (Figure 1c): It also develops hierarchically nested groundwater flow systems that include local, intermediate and regional flow systems under the control of the groundwater-fed lakes as discharge areas, as shown in Figure 2a [24,25]. The local flow system usually occurs between groundwater divides and neighboring lakes (e.g., Dake Lake and Mukai Lake) with a relatively shallow circulation depth, while, intermediate and regional flow systems are present between the major divides and the relatively larger lakes such as HTCH Lake, with significantly deeper circulation depths and longer flow distances relative to local flow systems. As shown in Figure 2a, the study area (HTCH Lake Basin) located in the lower reaches of the basin probably develops local, intermediate and regional flow systems (Figure 2). According to the potentiometric contours map of the shallow groundwater (Figure 1d), groundwater flows along the cross section A′–A″ from northwest to southeast, which agrees well with the topography in the study area (Figure 1d).

3. Materials and Methods

3.1. Sampling and Analysis

An important sampling campaign was conducted in the HTCH basin during September 2017. A total of 38 groundwater samples were collected along the cross section A′–A″ following the groundwater flow direction. Thirty-three samples were gathered from active water wells with depths from about 10 m to 200 m, and 5 samples were collected from deep boreholes with depths varying between about 350 m to 600 m. In addition, a total of 5 lake water samples were taken from the HTCH Lake. The sampling sites are shown in Figure 1c,d. The shallow groundwater samples are considered belonging to the same groundwater flow system, namely the local flow system, according to the hydrogeological context reported in previous study [14,15]. Meanwhile, in order to better understand the chemical evolution of groundwater, a total of 5 rock core samples with different depths were collected from a deep borehole, as shown in Figure 1c.

All groundwater samples were collected at a well head using pumps installed in the wells. In both cases, the wells were allowed to discharge for a period of time prior to sampling to ensure the removal of stored water from the well column. In situ cellulose membrane filters (0.45 μm, Jinteng, Tianjin Municipality, China) were used to filter the samples for analyses of water chemistry and Sr isotope; raw water samples without filtering were collected in sampling polyethylene bottles for hydrogen, oxygen, sulfur, carbon isotope and tritium determinations. The sampling bottles were prerinsed with water samples three times before the final water was gathered. Finally, all water sample bottles were sealed with parafilm tape in order to prevent them from evaporation.

Total dissolved solids (TDS), pH, electrical conductivity (EC), Eh and temperature parameters were measured in the field by a portable Multiparameter Digital Meter (HQ40D, Hach, Loveland, CO, USA) that was calibrated using standard solution before use. These field parameters were recorded in situ until they reached a steady state. The analyses for the groundwater performed in the laboratory included water chemistry and isotopes (δ²H, δ¹⁸O, ⁸⁷Sr/⁸⁶Sr, δ³⁴S, ¹⁴C and ³H). All lake water samples were measured for major ion chemistry (K⁺, Na⁺, Ca²⁺, Mg²⁺, Cl⁻, SO₄²⁻, HCO₃⁻, CO₃²⁻, and NO₃⁻). All groundwater samples were analyzed for major ion chemistry and stable water isotopes (δ²H and δ¹⁸O). Minor ion concentrations and Sr isotope ratios were determined for 15 representative samples. Ten samples were selected for the carbon isotopic analysis (δ¹³C and ¹⁴C), and tritium (³H) was measured for 7 representative samples. Major ions except HCO₃⁻ and CO₃²⁻ were analyzed using Ion Chromatography (Dionex ICS-2100, Sunnyvale, CA, USA) in the Water Isotope and Water–Rock Interaction Laboratory, Institute of Geology and Geophysics, Chinese Academy of Sciences (IGG, CAS). The alkalinity was measured by the titration method using 0.0048 M HCl immediately after returning the same laboratory. The charge balance error of major ions for all groundwater samples was within ±5%. The δ²H and δ¹⁸O compositions were analyzed using a liquid water isotope laser analyzer (Picarro L2120-i, Picarro, Inc., Santa Clara, CA, USA) in the same laboratory. The isotopic data was expressed in the standard delta notation (δ) as per mil (‰) relative to the Vienna Standard Mean Ocean Water (VSMOW, 0‰), and the measurement precision was ±0.5‰ for δ²H and ±0.1‰ for δ¹⁸O. The minor elements were measured by inductively coupled plasma mass spectrometry (ICP-MS NexION300D, PerkinElmer, Inc., Shelton, CT, USA) at the Beijing Research Institute of Uranium Geology, and the analytical precision is ±3% of the concentration based on the reproducibility of samples and standards. The ⁸⁷Sr/⁸⁶Sr ratios were determined using a thermal ionization mass spectrometer (Phoenix, ThermoFisher Scientific, Waltham, MA, USA) at the Beijing Research Institute of Uranium Geology with the analytical standard error (2σ) less than 0.00003. The δ³⁴S values for dissolved SO₄ in groundwater were measured by MAT 253 stable isotope mass spectrometry (ThermoFisher Scientific, Waltham, MA, USA) in the stable isotopes laboratory of IGG-CAS, and the analytical precision was approximately 0.2‰. The tritium (³H) was also measured at the Beijing Research Institute of Uranium Geology with the detection limit of 1.3 TU. The δ¹³C and ¹⁴C was determined in the Beta Analytic^TM in the United States by accelerator mass spectrometry (AMS) with an analytical precision of ±0.1‰ for δ¹³C and 0.1 pMC (percent Modern Carbon) for ¹⁴C. In addition, the aquifer mineralogical compositions were determined using an X-ray Diffraction (XRD) method at Lanzhou Center for Oil and Gas Resource, Institute of Geology and Geophysics, Chinese Academy of Sciences.

3.2. Inverse Geochemical Modeling

Geochemical modeling provides a useful tool for the quantitative evaluation of the chemical evolution in a particular groundwater system [26,27]. Inverse geochemical modeling is widely used to quantitatively simulate the key geochemical processes that can explain the differences observed between the initial and final water compositions along the flow path in a groundwater system [28,29]. In this study, the inverse geochemical modeling was conducted by means of the geochemical modeling program PHREEQC (version 2) developed by the United States Geological Survey (USGS) [30]. Furthermore, the saturation indices (SI) for minerals were also calculated with the help of this program.

4. Results and Discussion

The well locations, field data and major ion chemistry of groundwater and lake water are summarized in Table S1 in the Supplementary Material (SM). Table S2 in the SM presents the minor ion chemistry, isotopic compositions (δ²H, δ¹⁸O, ⁸⁷Sr/⁸⁶Sr, δ³⁴S, ³H, and ¹⁴C) and SI, along with the calculated ¹⁴C age of groundwater in the study area.

4.1. Hydrochemical Pattern Resulting from Groundwater Flow Systems

Based on the available chemical data, the hydrochemical characteristics of HTCH lake water and groundwater were discussed.

The chemical data for lake water show that TDS, a comprehensive indicator of hydrochemistry, is relatively high and ranges from 4464.5 mg/L to 5165.7 mg/L with the mean value of 4719.2 mg/L. The pH values lie between 9.97 and 10.03 with an average of 10.00, indicating an alkaline nature. The water type is Na–Cl(HCO₃), as shown in Piper plot (Figure 3). Na⁺ and Cl⁻ are the dominant elements: Na⁺ accounts for, on average, 88.6% of total cations, and Cl⁻ accounts for, on average, 51.3% of total anions. This is in accordance with hydrochemical characteristics of inland salt lake in other regions [16].

The TDS of groundwater range between 171.8 and 724.0 mg/L with mean and median values of 300.7 mg/L and 234.9 mg/L, respectively. Groundwater is slightly alkaline and has pH values varying from 7.69 and 9.39. The measured Eh values of groundwater range from −137.50 to −55.60 mV with the mean value of −90.34 mV. Ca²⁺, and the Na⁺ concentrations in groundwater vary in a wide range from 2.6 to 64.3 mg/L and 14.9 to 225.3 mg/L, respectively. The Cl⁻ and NO₃⁻ concentrations also exhibit a wide range from 3.5 to 60.9 mg/L and from ND (not detected) to 67.8 mg/L, respectively.

The Piper plot is widely used to graphically display the bulk major ion compositions of groundwater, and it is, therefore, used to classify groundwater [31]. Figure 3 shows that they are found to be chemically classified into five distinct types: Ca(Mg)–HCO₃ (type 1), mixed cations-HCO₃ (type 2), Na–HCO₃ (type 3), Na–HCO₃(SO₄) (type 4) and Na–SO₄ (type 5). In this classification scheme, type 1 through type 5 represent a natural chemical evolutionary sequence (Figure 3). The TDS and major ion concentrations for the 5 water types are also show illustrated in the Box plots of Figure 4. The calculated average values of the geochemical data for the 5 types of groundwater, including TDS, temperature, pH, Eh, major ions and stable isotopes (δD, δ¹⁸O and ⁸⁷Sr/⁸⁶Sr) are summarized in Table S3. The hydrochemical characteristics of each type of groundwater are discussed in detail below.

Type 1 waters have a Ca(Mg)–HCO₃ composition and average TDS of 234 mg/L (Table S3). These waters contain more than 70% (in meq/L) Ca²⁺ + Mg²⁺ of the total cations content. HCO₃⁻ is the major anion component, which accounts for more than 60% of total anions contents. This type predominantly occurs in the topographic highs of the basin, as shown in Figure 2b. Type 2 waters are mixed cations-HCO₃, and no dominant cation exists, belonging to the transitional type. The anion relationship is rHCO₃⁻ > rSO₄²⁻ > rCl⁻. Type 3 is Na–HCO₃ water. Na⁺ is the major cation, which constitutes more than 54% of the total cations; HCO₃⁻, the major anion, accounts for more than 67% of the total anions. These waters are low in SO₄²⁻ content that makes up less than 23% of the total anions. Type 4 is Na–HCO₃(SO₄) water with the mean TDS of 305.2 mg/L. Na⁺ is the predominant cation, and it accounts for greater than 77% of the total cation contents; HCO₃⁻ is the predominant anion, and it occurs in amounts greater than 48% of the total anion contents. The SO₄²⁻ constitutes more than 23% of the total anions, which is considerably higher relative to type 3 waters. These three types, including type 2, 3 and 4, are predominantly located in the topographic lows, as shown in Figure 2b. Type 5 is Na–SO₄ water, with the mean TDS value of 692.8 mg/L. Na⁺ contents range from 213.8 to 258.0 mg/L with the average value of 225.03 mg/L and constitute, on average, about 91% of total cations; SO₄²⁻ concentrations vary between 322.2 and 608.7 mg/L with the mean of 387.11 mg/L, which accounts for, on average, 68% of total anions. The concentrations of TDS, Na⁺ and SO₄²⁻ in type 5 are significantly higher than types 1, 2, 3 and 4, as illustrated in Figure 4. In summary, type 1 is characterized by and dominated by Ca²⁺, while types 3, 4 and 5 are characterized by Na⁺ among cations; types 1, 2, 3 and 4 are characterized by a predominance of HCO₃⁻, whereas type 5 is characterized by SO₄²⁻ among anions (Figure 3). In addition, it is should be noted that the mean Cl⁻ concentrations are much higher in type 1 (16.4 mg/L) than type 2 (9.0 mg/L) and type 3 (8.2 mg/L), as shown in Table S3 and Figure 4e, and is not consistent with the general pattern that Cl⁻ contents increase along the flow line. This is attributed to the fact that some samples in type 1 waters have been strongly affected by anthropogenic contamination, such as agricultural fertilizer, as indicated by higher mean NO₃ contents in type 1 (range: 13.6–67.8 mg/L; mean: 36.1 mg/L) than type 2 (range: 3.1–31.1 mg/L; mean: 18.5 mg/L) and type 3 (range: 13.6–180.1 mg/L; mean: 2.4 mg/L) (Table S3 and Figure 5). A similar contamination phenomenon caused by human activities is also observed in the regions close to the study area in the Ordos Plateau [20,32].

Broad relationships between the hydrochemical pattern and groundwater flow systems have been demonstrated by various workers. According to theory and observations, the groundwater flow from the recharge to discharge zones may be attended generally by an increase in the TDS and a systematic variation in the dominant cations and anions [33,34,35]. On the other hand, in a complex basin with the development of three types of groundwater flow systems (i.e., local, intermediate and regional), as described by Toth [36], groundwater from local flow systems may exhibit significantly distinct hydrochemical characteristics from those from intermediate/regional flow systems. This is because of the significant difference in the length of the flow paths and, hence, the contact time between groundwater and rock, especially in the lower reaches of a basin where a regional flow system lies below a local system [37,38]. Generally, systems may be characterized as local by low-TDS bicarbonate waters, as intermediate by intermediate-TDS-sulfate waters and as regional by high-TDS chloride-rich waters.

When we correlate the water types with the geographic locations, types 1, 2, 3 and 4 characterized by a predominance of HCO₃⁻ were found to be located in the shallow aquifer (<200 m), whereas all samples of type 5 dominated by SO₄²⁻ occurred in the deep parts (around 350–600 m), as shown in Figure 2b and Figure 5. Such a vertical pattern of groundwater chemistry is in accordance with the propositions of Freeze [37] and Toth [38]. In addition, groundwater is of a meteoric origin, and its chemical composition is predominantly controlled by a water–rock interaction within the aquifer in the study area that is relatively simple and macroscopically homogeneous sandstone (see Section 4.3). Therefore, the significant geochemical differences between shallow and deep groundwater seem to support the hypothesis of nested flow systems with different flow distances and flow times (Figure 2). That is, the groundwater chemistry suggests that shallow groundwater (types 1, 2, 3 and 4) belongs to a local flow system with a relatively short travel time and distance, while deep groundwater (type 5) belongs to a regional flow system with a longer travel time and distance. The TDS and Cl⁻ content, traditional indicators of groundwater residence time, are significantly higher in deep groundwater compared to shallow groundwater, as shown in Figure 5a,g, which also support this interpretation.

4.2. Groundwater Origin

The environmental isotopes of deuterium (²H) and oxygen-18 (¹⁸O) are excellent tracers for determining the groundwater origin [39]. A local meteoric water line (LMWL: δ²H = 6.45δ¹⁸O − 6.51) for the Ordos Plateau basin defined by Yin et al. [23] provides the basis for the interpretation of δ²H and δ¹⁸O values for groundwater in this study. The lower slope of LMWL relative to the global meteoric water line (GMWL) (δ²H = 8δ¹⁸O − 10) is attributed to the secondary evaporation process during rainfall as a result of an arid environment.

The δ²H and δ¹⁸O values in shallow groundwater exhibit a relatively wide range from −86.5‰ to −57.5‰ and from −11.3‰ to −7.4‰, respectively, which is within the general range previously observed for shallow groundwater in the Ordos Plateau [23,40]. The plot of δ²H versus δ¹⁸O (Figure 6) shows that all groundwater samples lie close to the LMWL and GMWL, indicating that they are of local meteoric origin. However, they can be divided into two groups. One group (group 1) is the type 1 groundwater located in the recharge area of local flow system; the average stable isotopic composition (−63.3‰ for δ²H and −8.3‰ for δ¹⁸O) in these groundwater corresponds to that of modern rainfall in the Ordos Plateau with a weighted mean of −58.5‰ for δ²H and −8.1‰ for δ¹⁸O (Figure 7) [23]. This suggests that they are recharged primarily by modern precipitation or precipitation in the late Holocene period. The other group (group 2) includes the other four types of groundwater (types 2, 3, 4 and 5). They are depleted in heavy isotopes (mean −80.9‰ for δ²H and −10.4‰ for δ¹⁸O) in comparison with type 1 groundwater and present meteoric water (Figure 7), suggesting that they are primarily recharged during the late Pleistocene and early Holocene periods when climate was relatively colder and wetter than the present, after having considered that the altitude effect was not significant in the whole Ordos Plateau [23,40]. In addition, in group 2, type 5 waters are found to have isotopic compositions very similar to type 4 and eight samples in type 3 because they are both old waters, while type 2 waters and four samples in type 3 display the isotopic compositions within those of type 1 and those 8 samples in type 3 (type 4, 5), which may be caused by a mixing between type 1 and type 3 based on the chemical characteristics and hydrogeological setting.

4.3. Mechanisms Controlling Groundwater Chemical Evolution

This section aims to identify the dominant mechanisms responsible for the chemical compositions of groundwater of various types. Understanding geochemical processes in groundwater is a prerequisite for the correction of ¹⁴C ages and for recording groundwater recharge and flow. Processes controlling the hydrochemical evolution were determined using hydrochemistry, stable isotope and inverse geochemical modeling.

4.3.1. Constraints from Water Chemistry

Generally, the chemical compositions of groundwater are primarily controlled by various factors, including the chemical composition of recharge water, water–rock interaction, dissolution/precipitation of mineral phases, evaporation, mixing processes and human activities [1,2]. The Gibbs diagram [41], which exhibits the ratios of Na²⁺/(Na²⁺ + Ca²⁺) and Cl⁻/(Cl⁻ + HCO₃⁻) as a function of TDS, is an effective tool to assess the natural mechanisms controlling water chemistry, including the atmospheric precipitation, the water–rock interaction and the evaporation-crystallization process. Figure 8 shows that all groundwater samples fall into the zone of rock dominance, suggesting that water–rock interactions are the major mechanism controlling groundwater chemistry. The Cretaceous sandstones consist mainly of quartz; K-feldspar; plagioclase; calcite; and minor amounts of dolomite, ankerite, siderite, pyrite, hematite, analcite, Glauber’s salt and clay minerals, as presented in the Table S4. Therefore, potential water–rock interactions include the dissolution/weathering of these minerals forming aquifer and cation exchange occurring on the clay minerals. In addition, the dissolution of halite and gypsum can also take place in the aquifer, although they are not found in the study area. This is due to the fact that they have been observed in the Cretaceous aquifer in the Ordos Plateau [15].

In order to further reveal the origin of solutes in groundwater, the relationship between the dissolved chemical constituents is employed [42]. The relationship between Na⁺ and Cl⁻ is commonly used as a first order indicator of water–rock interactions [43]. The plot of Na⁺ versus Cl⁻ (Figure 9a) shows that only three samples in type 1 groundwater samples (about 25% of type 1 waters) lie close to the 1:1 relation line expected from the dissolution of halite, indicating the occurrence of halite dissolution in these waters, while other samples cluster above the 1:1 line and molar Na/Cl ratios vary between 1.5 and 30.2 with a mean value of 8.7. This indicates the existence of additional sources contributing to the dissolved Na⁺, which potentially include the weathering of Na-bearing silicates, a dissolution of Na-bearing evaporites and cation exchange between Ca²⁺/Mg²⁺ in groundwater and Na⁺ on clay minerals after taking into consideration the mineralogical compositions of aquifer studied [44,45,46].

The plot of (Ca²⁺ + Mg²⁺) versus (HCO₃^– + SO₄^2–) (Figure 9b) indicates that type 1 groundwater is distributed around the 1:1 relation line, suggesting that Ca, Mg and HCO₃ in these waters are derived mainly from simple dissolution of calcite, dolomite and gypsum. The mean SI values for calcite and dolomite are 0.71 and 0.93, respectively, indicative of a slight oversaturation of groundwater with carbonates minerals and the importance of carbonates dissolution in these waters. These reactions commonly occur in the recharge areas in many sedimentary aquifer environments, and the reaction equations are as follows [2]:

CaCO₃ + CO₂ + H₂O = Ca²⁺ + 2HCO₃⁻

(1)

CaMg(CO₃)₂ + 2CO₂ + 2H₂O = Ca²⁺ + Mg²⁺ + 4HCO₃⁻

(2)

CaSO₄ = Ca²⁺ + SO₄²⁻

(3)

In contrast, type 2, 3, 4 and 5 waters exhibit a significant deviation from the 1:1 relation line and plot above it, indicating the simple dissolution of carbonate and gypsum alone could not interpret their chemical compositions. This may be caused by the cation exchange, silicates weathering or evaporites dissolution [47]. Figure 9c shows that majority of these waters plot above or below the theoretical line of weathering of plagioclase (y = x), suggesting silicate weathering is not the major source of Na⁺. As shown in Figure 9d, type 2, 3, 4 and 5 waters have a much higher Na/Ca ratio relative to type 1 water and TDS slightly increases as groundwater evolves from type 1 to types 2, 3 and 4. This seems to imply that a cation exchange has taken place as groundwater flows in the aquifer. The plot of (Na⁺ − Cl⁻) against (Ca²⁺ + Mg²⁺) could be employed to evaluate the cation exchange [46,48]. If the relationship between these two parameters is linear with a slope of −1.0, suggesting a cation exchange is significant in controlling geochemical compositions. As shown in Figure 9e, types 1, 2, 3 and 4 are located near the 1:1 relation line; in contrast, type 5 groundwater lie away from the 1:1 line. This mirror image of mono- (Na⁺) and bivalent cations (Ca²⁺, Mg²⁺), along with the almost constant HCO₃^– contents along the groundwater flow path (Figure 5e), indicate that a cation exchange controls the hydrochemical evolution of shallow groundwater (types 1, 2, 3 and 4) along the flow path; the reaction equations are as follows:

NaX + Ca²⁺/Mg²⁺ = CaX₂/MgX₂ + Na⁺

(4)

The use of two Chloro-Alkaline indices (CAI-1 and CAI-2) proposed by Schoeller is helpful for indicating the cation exchange process [49]. CAI-1 and CAI-2 are defined as the ratio of (Cl − (Na + K)) to Cl and (HCO₃ + SO₄ + CO₃ + NO₃), respectively. The positive indices suggest that the ion exchange occurs between Na⁺ in the groundwater and Ca²⁺/Mg²⁺ in the aquifer sediments, while the negative values suggest that the reverse ion exchange takes place. As shown in Figure 9f, both of the indices (CAI-1 and CAI-2) are found to be below zero in the type 2, 3 and 4 groundwater samples, implying that a high concentration of Na⁺ and a low concentration of Ca²⁺ are attributed to the ion-exchange process in these waters.

Figure 9g shows that deep type 5 water displays a significant deviation from the 1:1 line expected from gypsum dissolution and plots above the 1:1 line, indicating that gypsum is not the major process in this type of water and that SO₄²⁻ has additional sources. The SI values for gypsum in these waters range from −1.85 to −1.67 with an average of −1.79, furthering supporting this interpretation. However, these waters have a relatively higher average SI value for gypsum than type 1 (−2.33), type 2 (−2.46), type 3 (−2.84) and type 4 (−2.84), as shown in Table S3, suggesting a stronger degree of gypsum dissolution in this type than in the other four types. Figure 9h shows that type 5 groundwater plots relatively close to the 1:1 relation line, suggesting Na⁺, SO₄²⁻ and dominant ions in type 5 primarily derive from the dissolution of evaporite including halite and Glauber’s salt (Na₂SO₄). The reaction equation of Glauber’s salt dissolution is as follows:

Na₂SO₄ = 2Na⁺ + SO₄²⁻

(5)

Another significant characteristic of the groundwater is the extremely low K⁺ content between 0.6 and 4.4 mg/L (mean 2.3) and the very low K/Na molar ratios of near 0. The low levels of potassium may be attributed to its tendency to participate in the formation of secondary minerals and the lower reactivity of K-feldspars compared to plagioclase [8,43].

The dissolved silica concentrations (range: 1.66–11.90 mg/L; mean: 7.14 mg/L) are generally lower than the typical case where silicate weathering is a major process in groundwater (approx. 20 to >100 mg/L) [9], suggesting that the silicate weathering is a minor source of Na⁺ in groundwater.

4.3.2. Constraints from Strontium Isotope

The Sr isotopes are extensively applied as tracers to effectively delineated mineral weathering reactions in subsurface hydrology [50,51,52]. This is attributed to the following properties. The strontium can be removed from the water through mineral precipitation or cation exchange reaction, but the Sr isotopic composition of water remains largely unaffected because its mass-dependent fraction is usually negligible in nature. Moreover, Sr²⁺ is geochemically similar to Ca²⁺ and usually resides in carbonates, gypsums and K-bearing or Ca-bearing silicates [51,53]. These common rock-forming minerals generally exhibit a wide and predictable range of ⁸⁷Sr/⁸⁶Sr ratios that are controlled by initial values, Rb/Sr ratios and the age of minerals. Silicate minerals, particularly potassic silicates (e.g., biotite and K-feldspar), that exceed a few million years old have more radiogenic ⁸⁷Sr/⁸⁶Sr ratios due to the addition of radiogenic ⁸⁷Sr through the decay of ⁸⁷Rb (t_1/2 = 4.8 × 10¹⁰ years). In contrast, Ca-rich minerals such as calcite or gypsum display relatively less radiogenic ⁸⁷Sr/⁸⁶Sr ratios and remain relatively constant over time.

Groundwater primarily originates from the local precipitation in the study area, as indicated by δD and δ¹⁸O. Therefore, the strontium geochemistry of groundwater is initially controlled by the precipitation. The Sr concentrations in precipitation samples collected from the Otak meteorological station (Figure 1b), nearest to the study area, range from 10 to 250 μg/L with the volume-weighted mean (VWM) value of 50 μg/L [54]. The Sr concentrations in groundwater vary in a wide range from 377 to 1000 μg/L, with an average value of around 630 μg/L. The much higher Sr²⁺ concentration in groundwater compared to those observed for local precipitation demonstrates that Sr²⁺ in groundwater have other potential sources besides precipitation, which may include the dissolution of carbonate and silicate minerals that are major Sr-bearing minerals in the aquifer matrix.

The Sr sources might control the ⁸⁷Sr/⁸⁶Sr ratios and can be well constrained by the relationship between the Sr isotopic compositions and concentration ratios for groundwater [55]. Figure 10a shows their relationship in the plot of ⁸⁷Sr/⁸⁶Sr ratios versus 1/Sr. The geochemical compositions of three potential end-members are also shown: (a) carbonate fraction of aquifer rocks, represented by an ⁸⁷Sr/⁸⁶Sr ratio and Sr concentrations of 0.710148 and 309 μg/g, respectively [56]; (b) silicate fractions of aquifer rocks, characterized by a significantly more radiogenic ⁸⁷Sr/⁸⁶Sr ratio (0.716552) and a relatively lower Sr content (284 μg/g) relative to carbonate minerals [57]; and (c) local precipitation, characterized by an ⁸⁷Sr/⁸⁶Sr ratio and Sr concentration of 0.710760 and 309 μg/g, respectively [54]. Figure 10a shows almost all the groundwater samples points are located within the triangle consisting of local precipitation, carbonate and silicate fractions, and the plot is obviously closer to the carbonate fraction end-member than the other two end-members, which indicates the carbonate dissolution exerts the significantly more important control on the Sr²⁺ in groundwater than the other two sources. Figure 10b shows that ⁸⁷Sr/⁸⁶Sr ratios tend to gradually decrease with increasing Sr²⁺ concentrations in groundwater, furthering demonstrating that the major source of the dissolved Sr²⁺ is not the weathering of silicates with more radiogenic ⁸⁷Sr/⁸⁶Sr ratios. This is due to the fact that ⁸⁷Sr/⁸⁶Sr ratios will increase with increasing Sr²⁺ concentrations if silicate weathering is an important process controlling the groundwater hydrochemical evolution [9]. Figure 10b displays that types 3, 4 and 5 have relatively lower ⁸⁷Sr/⁸⁶Sr ratios compared to type 1, which may result from the gypsum dissolution as groundwater moves from recharge areas within the aquifer. This is supported by the fact that the mean SO₄²⁻ concentrations are higher for type 3 (28.6 mg/L), type 4 groundwater (67.1 mg/L) and type 5 (342.8 mg/L) than type 1 groundwater (18.9 mg/L) as shown in Table S3.

As the silicate weathering is not significant within the aquifer, excess Na⁺ relative to the Cl⁻ concentration in type 3 and 4 groundwater is believed to be the primary result of the cation exchange process and the dissolution of Glauber’s salt. Types 3 and 4 are characterized by a predominance of Na⁺ and HCO₃⁻; therefore, cation exchange is the major source of excess Na in these waters. That is, Sr geochemical signatures further confirm the conclusion that a cation exchange significantly contributes to dissolved Na⁺ in type 3 and type 4 groundwater, deduced from the groundwater chemistry discussed above.

4.3.3. Constraints from Sulfur Isotope

Groundwater sulfate can be derived from various sources, which include the accession of inorganic SO₄ via the atmospheric deposition of marine SO₄ aerosols, the dissolution of SO₄ minerals such as gypsum, the oxidation of reduced sulfide minerals such as pyrite or the mineralization of organic S compounds within the soil zone [58,59]. The sulfur isotopes values for SO₄ (δ³⁴S_SO4) of these sources are markedly different. The δ³⁴S values of sedimentary reduced sulfides (e.g., pyrite) typically range from approximately −50 to +10‰, while the evaporite minerals tend to be relatively enriched in δ³⁴S values between around 10 and 20‰ [13]. Therefore, the sources of SO₄²⁻ in groundwater can be better determined by sulfur isotope values of SO₄²⁻ [60,61].

In this study, an isotopic (δ³⁴S_SO4) analysis was carried out for type 5 groundwater (Na–SO₄) with relatively higher SO₄²⁻ concentrations in order to identify the SO₄²⁻ source in groundwater. Groundwater primarily originates from the local precipitation in the study area. Hence, the sulfur geochemistry of groundwater is initially controlled by the precipitation. Precipitation has SO₄²⁻ concentrations between 9.56 and 28.67 mg/L, and the mean δ³⁴S value for SO₄ in precipitation is around 3.7‰ in the Ordos Plateau [62]. The SO₄²⁻ concentrations range from 322.2 to 367.4 mg/L in type 5 groundwater, which is significantly higher compared to local precipitation. This indicates the presence of additional SO₄²⁻ sources in groundwater except precipitation. The relatively higher δ³⁴S values of these waters (12.4‰ to 14.1‰, mean 13.0‰) than those of precipitation also support this conclusion. These waters have δ³⁴S values within the typical range of evaporite minerals, indicating that SO₄²⁻ in type 5 is primarily from the dissolution of evaporite such as gypsum and Glauber’s salt in the aquifer. This interpretation is consistent with the hydrochemical data discussed above.

4.3.4. Constraints from Inverse Geochemical Modeling

An inverse geochemical modeling was performed in order to further constrain and quantify the main natural geochemical processes identified above: the dissolution of calcite, dolomite, anhydrite, halite and Glauber’s salt; the precipitation of calcite and dolomite; and a cation exchange with clays. The weathering of silicate minerals is thought to be unimportant for the chemical evolution and, hence, is not considered in the models.

During the inverse modeling, at least two chemical analyses of groundwater at different points of the flow path and a set of plausible phases (mineral and/or gases) that potentially react along this path are required. In the models, the average chemical parameter values from the statistical analyses of each type of groundwater, as presented in Table S3, were used to represent the “initial” and “final” waters along the groundwater flow paths. A total of five models (models 1, 2, 3, 4 and 5) have been run for an interpretation of the chemical compositions of the five types of groundwater present in the study area. The five simulated paths are represented by the transition between the following end-members: (1) pure water to type 1 water; (2) type 1 water to type 2 water; (3) type 2 to type 3 water; (4) type 3 to type 4 water; and (5) type 4 to type 5 water. The phases used were calcite, dolomite, gypsum and halite in all the models based on the lithological and mineralogical information. CO₂ gas was only included in Model 1 from the precipitation to type 1 water to simulate carbonate dissolution under the action of CO₂ because only this system is supposed to be open for CO₂ in the recharge area. CaX, MgX and NaX phases were only included in Models 2, 3 and 4 to simulate the cation exchange; Glauber’s salt was only considered in Model 5 to simulate Glauber’s salt dissolution. The constraints used included Ca²⁺, Mg²⁺, Na⁺, HCO₃⁻, SO₄²⁻ and Cl⁻ in all models.

Table S5 presents the selected results of the inverse modeling from all possible models based on the statistical parameters (the sum of residuals and maximum fractional error) calculated by PHREEQC, representing different possible combinations of reactants and products that can account for the changes in the groundwater chemistry. This result demonstrates that the amount of a given mineral phase that must dissolve or the precipitation to produce the observed variations in groundwater chemistry between two end points along the flow path.

Model 1 is the geochemical path from the local precipitation to a typical recharge water of type 1 (Ca·Mg–HCO₃); the pure water is assumed to be representative of precipitation in this model. It is summarized as follows:

Precipitation (pure water) + CO₂ gas + Calcite + Dolomite + Gypsum + Halite→Ca·Mg–HCO₃ type water

Model 2 is the geochemical path from type 1 to type 2 (mixed cations-HCO₃) and is summarized as follows:

Ca·Mg–HCO₃ type water + Na from cation exchange + Gypsum→mixed cations-HCO₃ type water + calcite precipitation (or Ca loss for ion exchange) + dolomite precipitation (or Mg loss for ion exchange)

Model 3 reproduces the geochemical path from type 2 groundwater to type 3 (Na–HCO₃) and is summarized as follows:

Mixed cations-HCO₃ type water + Na from cation exchange + Gypsum→Na–HCO₃ type water + Ca loss for ion exchange + Mg loss for ion exchange

Model 4 reproduces the geochemical path from type 3 water to type 4 (Na–HCO₃·SO₄) and is summarized as follows:

Na–HCO₃ type water + Na from cation exchange + Gypsum + Halite→Na–HCO₃·SO₄ type water + calcite precipitation (Ca loss for ion exchange) + dolomite precipitation (Mg loss for ion exchange)

Model 5 represents the geochemical path from type 4 water to type 5 (Na–SO₄), and it is summarized as follows:

Na–HCO₃·SO₄ type water + Glauber’s salt + Gypsum + (Halite)→Na–SO₄ type water + calcite dolomite precipitation

Summarizing, the geochemical modeling results support the hypothetical geochemical processes derived from the qualitative analyses for the groundwater chemical and isotopic compositions. From the modeling results, it can be confirmed that the dissolution of carbonate and evaporite minerals, as well as a cation exchange within the aquifer primarily control the chemical evolution of groundwater in the study area.

4.4. Groundwater Age

Groundwater age can be used to further constrain groundwater origin and flow systems. Radiocarbon (¹⁴C) is the most widely applied dating method by far because of the almost ubiquitous presence of dissolved inorganic carbon (DIC) in groundwater [13].

In spite of the widespread use of ¹⁴C, the interpretation of ¹⁴C model age still is limited by many uncertainties in determining the ¹⁴C content of dissolved carbon in the recharge areas and in accounting for various chemical and physical processes affecting the ¹⁴C content along groundwater flow paths in aquifers [63,64,65,66]. Using the measurements of ¹⁴C of the DIC in groundwater, ¹⁴C_DIC, the basic equation for groundwater dating by using ¹⁴C is as follows [63]:

$$\mathrm{t}=-\frac{5730}{\ln 2}\ln (\frac{{}^{14}{\mathrm{C}}_{\mathrm{DIC}}}{{}^{14}{\mathrm{C}}_0})$$

(6)

where t is the groundwater age; 5730 is the modern half-life of ¹⁴C; ¹⁴C_DIC is the measured ¹⁴C value of the DIC in groundwater sample; and ¹⁴C₀ is the initial ¹⁴C content of DIC after correction for geochemical reactions without radioactivity decay.

The estimation of ¹⁴C₀ is essential for groundwater ¹⁴C dating. Numerous single-sample-based traditional geochemical adjustment models have been proposed to calculate ¹⁴C₀ so far, such as Vogel, Tamer, Pearson, Fontes and Garnier models and so on [67,68,69,70]. These different models may give significantly distinct values of ¹⁴C₀. The incorrect use of the models for chemical and isotopic measurements in groundwater will lead to a wide range of estimations of ¹⁴C₀ and greatly limits the usefulness of radiocarbon as a dating tool for groundwater [63]. Therefore, choosing the most appropriate models is of critical importance for improving the estimation of ¹⁴C age in a particular groundwater system, depending on which models most completely account for the geochemical processes affecting the chemical and isotopic compositions of the waters [65].

The Han and Plummer diagram of δ¹³C versus ¹⁴C [63] can aid in determining which adjustment model is the most appropriate and should be applied for the specific groundwater system. Figure 11 shows a plot of ¹⁴C versus δ¹³C (Han and Plummer diagram) for the groundwater samples studied, indicating that all samples can be divided into two groups. One group is composed of type 1 waters; they plot close to the IAEA (International Atomic Energy Agency) line, suggesting they are young waters with ¹⁴C age close to zero. Another group includes type 4 and 5 waters, which have δ¹³C_DIC ranging between δ¹³C_i (approx. −11.05‰) and δ¹³C_E (approx. −0.15‰) and are plotted in the Eichinger area. This indicates that the carbon isotopic compositions of DIC in these waters have been affected by a carbon exchange between DIC in groundwater and solid carbonate under closed system, and the Eichinger model [71], accounting for a carbon exchange between solid carbonate and DIC under a closed system, should be applied to estimate ¹⁴C₀ for type 4 and 5 waters. In fact, the models of Wigley and Evans also take into consideration the carbon exchange between DIC and carbonate mineral, but they are less precise compared to the Eichinger model because they neglect the contribution of CO₂ (aq) to DIC (the models assume that DIC = HCO₃⁻). Therefore, Eichinger’s model considering the contribution of CO₂ (aq) to DIC is selected to obtain a relatively accurate age in this study.

According to Han and Plummer, the Han and Plummer’s model and IAEA model are also applicable for groundwater samples studied [63]. Han and Plummer’s model is a revised version of the Fontes and Garnier model [72], which considers not only carbon mixing but also carbon exchange occurring both in the saturated and unsaturated zone. This model is a combination of three models: Tamers’, Mook’s and Eichinger’s models. If the carbon exchange between DIC and solid carbonate occurs but without an exchange with soil CO₂, the models reduces to Eichinger’s model. The carbon exchange between DIC and solid carbonate predominantly under closed system conditions occurs in the groundwater system studied; therefore, the Han and Plummer model will give a similar result to Eichinger’s model. Unlike the Han and Plummer’s model, assuming either open or closed systems, the IAEA model assumes mixed open and closed system evolutions [73]. That is, this model assumes that the isotopic composition of DIC evolved initially under completely open system condition and then evolved under closed system condition.

The measured and assumed parameters for estimating ¹⁴C₀ using various adjustment models are summarized in Table S6. The measured parameters include temperature, pH, DIC, contents of dissolved CO₂ (aq), HCO₃⁻ and CO₃⁻. The assumed parameters include ¹⁴C and δ¹³C values of soil gas CO₂ and solid carbonate mineral. The ¹⁴C of soil CO₂ (¹⁴C_g) and carbonate mineral (¹⁴C_s) are assumed to be 100 pMC and 0 pMC, respectively. The δ¹³C value for carbonate mineral (δ¹³C_s) is assumed to be 0‰. The vegetation in the Ordos Plateau is predominantly C3 species, and the soil profiles in Yulin, approximately 125 km away from study area, show that δ¹³C values range from −22‰ to −20‰ with an average of −21‰ [74,75]. Therefore, the δ¹³C value for soli gas CO₂ (δ¹³C_g) is assumed to be −21‰ in the study area. In addition, the calculated fraction factor for determination of ¹⁴C₀ is also presented in Table S6. The ε_a/g and ε_b/g are ¹³C fractionate factors of dissolved CO₂ and HCO₃⁻ with respect to gaseous CO₂, respectively. The ε_a/s and ε_b/s are the ¹³C fractionate factors of dissolved CO₂ and HCO₃⁻ with respect to carbonate mineral.

The measured ¹⁴C_DIC values ranged from 1.4 to 89.6 pMC, with an average of 34.9 pMC (Table S7). Based on above these parameters, the calculated ¹⁴C₀ and ¹⁴C age including the apparent and corrected ages are presented in Table S7. The models of Eichinger and of Han and Plummer yield similar ages, which are relatively lower than those determined by IAEA model. The results derived from these three models are overall comparable, suggesting that ¹⁴C method provides the relatively reliable and reasonable ages to represent the mean groundwater residence time. The average values of three models were considered as the final ages. For shallow groundwater in the recharge areas (type 1), the corrected ¹⁴C age is approximately 0 year BP (Before Present); for those in discharge areas (type 3 and 4), the corrected ¹⁴C age ranges from 3929 to 16,996 years BP with the mean value of 9450 years, which is consistent with the ³H content that is below the detection limit of 1.3 TU (Table S2). For deep groundwater (type 5), the corrected ¹⁴C age varies from 23,145 to 27,436 years BP with the average of 25,075 years, which is also in agreement with ³H content below the detection limit (Table S2). Summarizing, the age data further demonstrates that shallow groundwater in the recharge areas were recharged during the Holocene period, while those in the discharge areas and deep groundwater were recharged during the late Pleistocene period. That is, groundwater age agrees well with the δD and δ¹⁸O data discussed above. In addition, it is found that deep groundwater is about 2.7 times higher than shallow groundwater in the discharge areas. It is well-known that the groundwater age should increase along the groundwater flow direction. Therefore, such a vertical distribution of groundwater age further corroborates the hypothesis on the development of hierarchically nested groundwater flow systems. That is, groundwater age indicates that shallow and deep groundwater belong to local and regional flow system, respectively.

It is well-known that ¹⁴C ages are not identical to calendar age because atmospheric ¹⁴C concentration changes temporarily as a result of variation in the production rate, resulting from geomagnetic, solar modulation of the cosmic-ray flux and the carbon cycle [63,76]. In order to convert ¹⁴C ages into calendar ages relative to 1950 AD, many radiocarbon age calibration curves have been constructed using a variety of methods. In this study, the calibration curve IntCal13, representing the midlatitude Northern Hemisphere atmospheric reservoirs, was applied to estimate calendar age. Before it is applied, a conversion of conventional ¹⁴C ages calculated by using the modern ¹⁴C half-life (5730 years) to the Libby half-life (5568 years) is required because the calibration curves are obtained by first using the Libby half-life. The calibration curve IntCal13 and groundwater samples are shown in Figure 12, and the average calendar years of groundwater are also summarized in Table S7. The results of the calendar ages for groundwater show that shallow groundwater in the recharge areas is approximately 0 year BP, while those in discharge area ranges from 4.340 to 19,965 years BP with the mean value of 10,829 years; deep groundwater varies from 26,736 to 28,125 years BP with the average of 28,588 years.

5. Conclusions

Groundwater is the main water resource in the Ordos Plateau, which is located in arid-semiarid regions of China. In this study, environmental tracers including hydrochemistry and isotope, along with the inverse geochemical modeling technique were integrated successfully to obtain a comprehensive understanding on the origin and hydrochemical evolution of shallow and deep groundwater in the HTCH Lake Basin, NE Ordos Plateau. The major conclusions are summarized as follows:

1. The geochemical data support the hypothesis of hierarchically nested groundwater flow system and suggest the development of local and regional flow systems in the study area. Based on the isotopic (δ²H, δ¹⁸O, ³H and ¹⁴C) compositions, shallow and deep groundwater both originate primarily from local precipitation. Shallow local groundwater in the recharge areas and discharge areas was recharged during the Holocene and late Pleistocene periods, respectively; deep regional groundwater was recharged by precipitation in the late Pleistocene period. Therefore, from the viewpoint of sustainable use, shallow discharge zone groundwater and deep groundwater should not be exploited heavily.

2. The groundwater chemical pattern is closely related to the groundwater flow systems, and the hydrochemical evolution is primarily controlled by geochemical processes. The shallow local groundwater chemistry systematically changes along the flow path: The water type progressively evolves from type 1 (Ca·Mg–HCO₃) in the recharge areas through type 2 (mixed cations-HCO₃) to type 3 (Na–HCO₃) and type 4 (Na–HCO₃·SO₄) in the discharge areas. This is primarily controlled by carbonates dissolution and cation exchange, which principally occur in the recharge areas and discharge areas of local flow system, respectively. In addition, human activities (e.g., agricultural fertilizer) also affect the chemical compositions of shallow recharge zone groundwater, particularly Cl⁻ and NO₃⁻ contents. Deep regional groundwater belongs to type 5 (Na–SO₄), which is dominated by evaporite dissolution (e.g., gypsum and Glauber’s salt). Although the aquifer studied consists mainly of silicate-dominated Cretaceous sandstone, Sr isotope provides the compelling evidence that silicate weathering is not important in the study area.