Groundwater origin, flow regime and geochemical evolution in arid endorheic watersheds: a case study from the Qaidam Basin, northwestern China

Groundwater origin, flow and geochemical evolution in the Golmud River watershed of the Qaidam Basin was assessed using hydrogeochemical, isotopic and numerical approaches. The stable isotopic results show groundwater in the basin originates from precipitation and meltwater in the mountainous areas of the Tibetan Plateau. Modern water was found in the alluvial fan and shallow aquifers of the loess plain. Deep confined groundwater was recharged by paleowater during the late Pleistocene and Holocene under a cold climate. Groundwater in the low-lying depression of the central basin is composed of paleobrines migrated from the western part of the basin due to tectonic uplift in the geological past. Groundwater chemistry is controlled by mineral dissolution (halite, gypsum, anhydrite, mirabilite), silicate weathering, cation exchange, evaporation and mineral precipitation (halite, gypsum, anhydrite, aragonite, calcite, dolomite) and varies from fresh to brine with the water types evolving from HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na to Cl-Na, Cl-K-Na and Cl-Mg type along the flow path. Groundwater flow patterns are closely related to stratigraphic control and lithological distribution. Three hierarchical groundwater flow systems, namely local, intermediate and regional, were identified using numerical modeling. The quantity of water discharge from these three systems accounts for approximately 83 %, 14 % and 3 %, respectively, of the total groundwater quantity of the watershed. This study can enhance the understanding of groundwater origin, circulation and evolution in the Qaidam Basin as well as other arid endorheic watersheds in northwestern China and elsewhere worldwide.

1 Introduction

Closed basins in arid and semiarid areas (e.g., the Great Artesian Basin and Murray Basin in Australia, Minqin Basin and Qaidam Basin in China, Death Valley in United States) have been the focus of attention due to their water scarcity, fragile ecology and rich mineral resources related to salt lakes (Edmunds et al., 2006; Lowenstein and Risacher, 2009; Love et al., 2013, 2017; Shand et al., 2013; Stone and Edmunds, 2014; He et al., 2015; Cartwright et al., 2017; Priestley et al., 2017a; Xiao et al., 2017). Groundwater plays a vital role in water supply, ecology maintenance, transportation of chemical components, and the formation of oil, gas reservoirs and mineral resources in these basins (Toth, 1980; Jiang et al., 2014; Jiao et al., 2015; Xiao et al., 2017). Understanding the regimes of groundwater recharge, flow and hydrogeochemical evolution is essential to maintain proper management and implement sustainable utilization of groundwater and mineral resources, as well as maintain the ecological environment (Cartwright et al., 2010a; Herrera et al., 2016).

Figure 1 Location of the study area (a) within China, (b) within the Qaidam Basin, and (c) details of sampling location and groundwater and physiographic zones within the study area. Please note that the above figure contains disputed territories.

In the arid northwest of China there are many closed basins such as the Tarim, Qaidam, Junggar and Minqin basins, in which the low-lying discharge areas are occupied by saline lakes, salt playas and salt crusts. The Qaidam Basin (Fig. 1a, b), the largest closed basin of the Tibetan Plateau, has the most plentiful number of salt lakes and salt playas and almost all varieties of salt deposits (Zheng et al., 1993), as well as rich oil and gas reservoirs (Tan et al., 2011; Ye et al., 2014). Considerable research has been conducted to provide support for water supply and mineral resource exploitation in the basin (Chen and Bowler, 1986; Vengosh et al., 1995; Lowenstein and Risacher, 2009; Li et al., 2010; Tan et al., 2011; Hou et al., 2014; Ye et al., 2014; Chen et al., 2017). However, most of the previous studies focused on the groundwater in the piedmont areas (Wang and Ren, 1996; Wang et al., 2010; Zhang, 2013; Hou et al., 2014; Su et al., 2015; Xu et al., 2017), material source of salt lakes (Vengosh et al., 1995; Lowenstein and Risacher, 2009; Tan et al., 2011; Chen et al., 2015) and the evolution of salt lakes (Chen and Bowler, 1986; Chen et al., 2017). The systematic understanding of regional groundwater regimes is still inadequate. This would limit the comprehensive planning and management of groundwater and salt lake mineral resource exploitation, and finally make it difficult to safeguard the circulation of the groundwater system and maintain the eco-environmental balance. Therefore, several attempts have been made to understand the regional groundwater regimes (Tan et al., 2009; Gu et al., 2017; Xiao et al., 2017), but very little research reported the circulation and evolution of groundwater from the mountain pass area to the central terminal lake area due to the notable difficulties in moving through and access the swamps on the lacustrine plain. This would greatly limit the full understanding of the role of hydrogeological processes in the basin.

Hydrogeological survey efforts have been undertaken in the Golmud River watershed of the basin since 2015 and have developed a better understanding of regional hydrogeological conditions. The main objective of this study is to assess the regional hydrogeological regime of closed basins in the arid northwest of China, using the Golmud River watershed as a case study. To achieve this aim, a comprehensive approach using environmental isotopes (²H, ¹⁸O, ³H, ¹³C, ¹⁴C) and hydrochemistry coupled with numerical simulation was performed. Stable hydrogen and oxygen can provide valuable information on the origin and recharge environment of groundwater, and radioactive isotopes such as ³H and ¹⁴C record the residence time of groundwater (Cartwright et al., 2007; Awaleh et al., 2017; Huang et al., 2017). Hydrochemical composition has recorded the recharge water characteristics, hydrostratigraphic information, geochemical interaction and other processes along the groundwater flow path (Redwan and Moneim, 2015; Verma et al., 2016; Love et al., 2017) and thus can be used to track groundwater evolution. Numerical simulation of groundwater flow is an essential tool to synthesize hydrogeological information and reveal groundwater flow patterns (Bredehoeft and Konikow, 2012; Anderson et al., 2015; Tóth et al., 2016). The combination of these approaches is robust to reveal groundwater origin, flow regimes, renewability, hydrochemical evolution and interaquifer mixing, as well as surface water and groundwater interactions, etc., in basins with complex hydrogeology or sparse monitoring data, and has been successfully applied in many basins such as the Great Artesian Basin and Murray Basin in Australia, Michigan Basin in US, Minqin Basin and Ordos Plateau in China, Stampriet Basin in Africa (Edmunds et al., 2006; Banks et al., 2010; Love et al., 2013, 2017; Stone and Edmunds, 2014; Su et al., 2016; Cartwright and Morgenstern, 2017; Petts et al., 2017; Priestley et al., 2017b).

The specific aims of the present study are to: (1) identify the recharge source of groundwater, (2) assess the regional groundwater chemistry characteristics, (3) determine the controlling mechanisms of hydrogeochemistry, (4) delineate regional groundwater flow patterns and (5) ultimately establish systematic regional groundwater regimes from the mountain pass to the terminal lake in the typical Golmud watershed of Qaidam Basin. This study would provide insights into the origin, recharge environment, flow regime and geochemical evolution of regional groundwater in arid endorheic watersheds of Qaidam Basin and provide reference for other arid closed basins in northwestern China as well as similar endorheic watersheds worldwide.

2 Study area

The Qaidam Basin is a large closed basin located on the northeastern margin of the Tibetan plateau, surrounded by the Qilian Mountains to the north, the Kunlun Mountains to the south and the Altun Mountains to the west (Fig. 1b). The study area, Golmud River watershed (GRW), is located in the southern part of the Qaidam Basin hosting the second largest river, the Golmud River, running from the Kunlun Mountains in the south to the low-lying depression in the north central part of the Basin (Fig. 1c). The Qarhan salt lake is the largest salt lake in China, located at the northern margin of the GRW, adjacent to the Golmud River's terminal Lake Dabusun. The third largest city on the Tibetan plateau, Golmud City, is also located in the GRW.

The outcropping stratigraphy of the GRW ranges from Proterozoic to Quaternary in age. The Quaternary strata are found in the mountainous areas to the south. These strata have undergone magmatic activity, uplift and tectonic movements, as well as intense weathering, resulting in massive material sources of sediments to the basin. The Quaternary deposits have thicknesses ranging from hundreds of meters in the piedmont area to thousands of meters in the low-lying depression (basin center) (Zheng et al., 1993; Chen et al., 2017). Field surveys found that salt crusts are formed on the ground surface in locations near the terminal areas of streams. Core drilling records also show many salt-bearing deposits such as halite, calcium, sulfate and sodium sulfate were observed throughout the strata (Chen and Bowler, 1986). The regional Quaternary aquifers in the basin vary from single unconfined gravel and sand layers with hydraulic conductivity (K) greater than 50 m d⁻¹ in the alluvial fan to multilayers of silt and clay with hydraulic conductivity (K) ranging from 0.1 to 0.001 m d⁻¹ in the low-lying depression (basin center). Three continuous aquitards (clay layers) are found in the basin at depths of 60, 290 and 450 m (Fig. 8), which have significant influences on confining groundwater flow (Shao et al., 2017).

The climate in the GRW is extremely variable, both spatially and temporally. Precipitation in the Kunlun Mountains is more than 200 mm per year, but less than 50 mm in the basin, and also presents a gradual decreasing trend from the piedmont area to the low-lying central depression. The potential evaporation is extremely high (>2600 mm per year). This hyper-arid climate results in aquifers in the basin that do not obtain effective recharge from the local precipitation. Groundwater in the basin is mainly recharged by Golmud River seepage through the riverbed in the alluvial fan and bedrock lateral inflow at the southern mountain front, and flows from the alluvial fan in the south to the basin center in the north (Fig. 1c). Much of the groundwater overflows as springs at the front of the alluvial fan due to the fining of sediments in the aquifer's downdip. The depth to groundwater is less than 3 m in most areas from the front of alluvial fan to the basin center, resulting in significant potentially evaporate loss of groundwater. The regional groundwater finally discharges to the terminal lake and undergoes large evaporate loss (Shao et al., 2017).

Based on the terrain, sediments and hydrogeological condition, the study area can be divided into 5 zones. Zone 1 is the Kunlun mountainous area, and Zone 2 is the alluvial fan plain at the Kunlun piedmont. Zones 3, 4 and 5 occur on the loess plain, where Zone 3 is the main groundwater overflow (discharge) zone of the watershed, and Zone 5 is the terminal lake zone (low-lying depression of the watershed) with salt crusts and playa. Zone 4 is the transition zone (middle to lower stream area of the watershed) between Zone 3 and 5.

3 Materials and methods

3.1 Hydrochemical and isotopic sampling and analytical methods

A total of 228 water samples were collected from GRW in 2015 and 2016, including 180 groundwater samples and 48 surface water samples (42 river water and 6 lake water samples) (Fig. 1c). Groundwater samples were collected from both shallow phreatic aquifers and deep confined aquifers. Surface water samples were obtained along the Golmud River, as well as from Lake Qarhan, Lake Dabusun and other small lakes in the low-lying depression area (basin center). In addition, 1 snow (snowmelt water) sample, 8 precipitation samples and 90 brine water samples (groundwater) with hydrogen and oxygen stable isotope data were obtained from China's stable isotope geochemistry database (http://210.73.59.163/isogeochem/, last access: 31 May 2018). The location of the snow and precipitation samples are shown in Fig. 1c. The detailed locations of the 90 brine water samples are not known, but it is known that all of these samples were collected from the Qarhan salt playa and Bieletan salt playa (Fig. 1c).

For groundwater sampling, all wells and boreholes, except those that were artesian, were pumped for several well volumes to remove the stagnant water in the wells and boreholes and monitored until the electrical conductivity (EC) of the pumping water was stable. The sampling procedure followed is described in Huang et al. (2016) and Chen et al. (2011). Samples for major element analysis were collected in two 250 mL high-density polyethylene bottles after filtration using 0.45 µm filter membranes (Huang et al., 2016). Samples for tritium (³H) and stable isotope (²H, ¹⁸O) analysis were collected in 500 and 50 mL glass bottles, respectively, that were filled to overflowing after rinsing and were sealed tightly. ¹³C and ¹⁴C samples were collected by adding BaCl₂ and CO₂-free NaOH to 120 L groundwater at pH = 12, obtaining BaCO₃ for dissolved inorganic carbon (DIC) analysis (Chen et al., 2011). The method used eliminates contact with the atmosphere in order to avoid CO₂ atmospheric contamination.

Parameters such as the water temperature (T), pH and EC were measured in the field with an in situ multiparameter instrument (Multi 350i/SET, Munich, Germany), and redox potential (Eh) was also determined in situ using a portable ORP tester (CLEAN ORP30 Tester, California, US). Major chemistry and isotopes (²H, ¹⁸O, ³H, ¹³C and ¹⁴C) of the sampled water were analyzed at the Laboratory of Groundwater Sciences and Engineering in the Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences (Shijiazhuang, Hebei Province, China). Major cations (K⁺, Na⁺, Ca²⁺, Mg²⁺) were measured by inductively coupled plasma–mass spectrometry (Agilent 7500ce ICP-MS, Tokyo, Japan). Total dissolved solids (TDSs) and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ were determined by gravimetric analysis and acid-base titration, respectively. Cl⁻ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ were analyzed using spectrophotometry (PerkinElmer Lambda 35, Waltham, MA, USA). The ionic charge balance of all samples were within 5 % difference. δ¹⁸O, δ²H and δ¹³C were measured by isotope ratio mass spectrometry using a Finnigan MAT 253, and δ¹⁸O, δ²H were reported relative to the Vienna Standard Mean Ocean Water (VSMOW) standard, and δ¹³C was reported relative to the Vienna Pee Dee Belemnite (VPDB). The analytical errors are ±0.2 ‰ for δ¹⁸O, ±1.0 ‰ for δ²H and ±0.5 ‰ for δ¹³C. The tritium content was determined using electrolytic enrichment and a liquid scintillation technique (Chen et al., 2011) with the precision of ±0.3 TU. The activity of ¹⁴C was analyzed by liquid scintillation counting (1220 Quantulus), and expressed as a percentage of modern carbon (pMC) with the precision of ±0.3 % (Su et al., 2018).

3.2 Two-dimensional groundwater flow numerical simulation

It is assumed that the variation of density and viscosity of waters could be neglected for calculations involving most of the flow system (Zone 1–4). For simplicity, groundwater in the terminal lake zone (Zone 5) is also regarded as being mainly driven by gravity. Thus the equation governing variably saturated groundwater flow is as follows (Richards, 1931):

$${\displaystyle \frac{\partial }{\partial t}}\mathit{\varphi }S=\mathrm{div}\left[\mathrm{K}\mathrm{\nabla }\mathrm{h}\right],$$

where ϕ is the porosity, S is the liquid saturation, K is the hydraulic conductivity (m d⁻¹), and h is the hydraulic head (m). The TOUGH2 code (Transport Of Unsaturated Groundwater and Heat), which has quite robust simulation capabilities, is used to numerically solve this equation (Pruess et al., 1999).

The cross section parallel to the main direction of groundwater flow in GRW (Figs. 1 and 8) was selected for the two-dimensional flow simulation. This section starts at the mountain pass and ends at the terminal Lake Dabusun, with an approximate length of 100 km. Boundaries were specified according to the hydrogeology condition. The southern lateral boundary and top boundary in the alluvial fan are defined as given flux boundaries, and the bottom boundary of the section is regarded as a zero flux boundary. The springs and evaporation are set as mixed boundaries. The lake boundary in the basin center is specified as a given head boundary.

An irregular discretization was conducted vertically to capture the variation of the water table near the ground surface and also implement an efficient simulation. Cells are presented with the minimum thickness of 0.1 m near the ground surface and gradual increasing thickness downward, with a maximum thickness of about 80 m. Equal discretization was applied in the horizontal direction with a horizontal size of 1000 m for one cell. The initial permeability of various lithologies is specified based on the borehole drilling records and pumping test results, with the K_h (horizontal hydraulic conductivity) in the range of 10²–10⁻³ m d⁻¹ and anisotropy ratio $K_{\mathrm{h}}/K_{\mathrm{v}}=\mathrm{5}$–10 (K_v is vertical hydraulic conductivity) (Shao et al., 2017). In this study, the model is used to present the flow pattern under equilibrium conditions, and thus the recharge rates and hydraulic heads are given according to the annual average values. Evaporation was modeled using a newly developed method described by Hao et al. (2016), and the initial potential evaporation specified is 2600 mm per year. Springs are simulated using the DELV module in TOUGH2, and the productivity index (PI) specified in the DELV module is calculated using the following equation (Pruess et al., 1999):

$$\mathrm{PI}={\displaystyle \frac{\mathrm{2}\mathit{\pi }\left(k\mathrm{\Delta }z\right)}{\ln (\sqrt{A/\mathit{\pi }}/r)+s-\mathrm{1}/\mathrm{2}}},$$

where Δz is the layer thickness (m), A is the grid block area (m²), r is the spring radius (m), and s is the skin factor. Annual average observed hydraulic heads are used as natural constraints for the model calibration.

Table 1 Statistical summary of physical and chemical parameters of the surface water and groundwater in the Golmud Watershed, Qaidam Basin, China.

RW: river water; SGW: shallow phreatic groundwater; DGW: deep confined groundwater; FLW: relatively fresh lake water; SLW: salt lake water.

4 Results

4.1 Hydrochemistry of surface waters and groundwaters

The statistical summaries of chemical analysis results for surface water and groundwater are presented in Table 1. River waters (RW) from the mountain pass (Zone 2) to the low-lying depression (Zone 5) are slightly alkaline, with a range in pH values from 7.94 to 9.45. Fresh lake water (FLW) L1, which was sampled from the freshwater lake (relatively fresh compared to other salt lakes) recharged directly by river water in the low-lying depression (Zone 5), is also slightly alkaline, with a pH value of 8.98. Samples from the salt lakes such as Lake Qarhan and Lake Dabusun are slightly acidic, with values range from 6.03 to 6.28. Groundwater in the study area is neutral to slightly alkaline. The shallow phreatic groundwater (SGW) shows an evolving trend from slightly alkaline to slightly acidic along the flow path with the pH varying from 9.34 to 6.03. However, deep confined groundwater (DGW) samples are all slightly alkaline, with pH values ranging between 7.83 and 8.69. The redox potential of SGW is in the range of 123–162 mV from alluvial fan to middle to lower stream areas (Zone 2, 3 and 4), suggesting an oxidation condition. The Eh values of DGW vary from 153 to 40 mV along the flow path (Zone 3 to 4), indicating the redox condition gradually evolves from a state of oxidation to reduction (Fig. 3e).

Figure 2 Piper diagrams of sampled surface water and groundwater. (a) Surface waters; (b) shallow phreatic groundwaters; (c) deep confined groundwaters from the Qaidam Basin, China (Red dashed lines and arrows indicate the direction of evolutionary flow systems).

Figure 3 Isotopic data and Eh vs. distance from the mountain pass along the groundwater flow paths. (a) ¹⁸O vs. distance, (b) ²H vs. distance, (c) ¹⁴C vs. distance, (d) ³H vs. distance, (e) Eh vs. distance.

Surface water and groundwater present distinct major solute chemistry across the study area. As shown in Table 1, the concentration of ions in RW demonstrates an increase along river flow paths, with TDS values varying from 393 to 2319 mg L⁻¹. The TDS value of FLW (L1) is much higher than that of RW in the low-lying depression (Zone 5), with the TDS value of 10 937 mg L⁻¹, while the salt lake waters (SLWs) have extremely high TDS values ranging from 339 098 to 403 758 mg L⁻¹. The dominant ions of RW are ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and Na⁺ with the concentration range of 184–215 mg L⁻¹ for ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and 63–92 mg L⁻¹ for Na⁺, respectively, in the alluvial fan area (Zone 2). These dominant ions gradually evolve to Cl⁻ and Na⁺with the concentration range of 655–1 776 mg L⁻¹ for Cl⁻ and 438–996 mg L⁻¹ for Na⁺ in the low-lying depression (Zone 5). FLW (L1) has the same dominant ions with RW in the low-lying depression (Zone 5), but with a higher concentration of 5912 mg L⁻¹ for Cl⁻ and 2957 mg L⁻¹ for Na⁺. SLW is dominated by Cl⁻ and Mg²⁺ with the concentration range of 276 849 to 285 780 mg L⁻¹ for Cl⁻ and 99 500 to 100 240 mg L⁻¹ for Mg²⁺, respectively. Overall, the surface water types evolve from HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na type in the alluvial fan area (Zone 2) to Cl-Na, Cl-K-Na and Cl-Mg types in the low-lying central depression (Zone 5) (Fig. 2a).

Table 2 Statistical summary of isotopic analysis results of precipitation, surface water and groundwater in the Golmud Watershed, Qaidam Basin, China.

SNW: snowmelt water; PW: precipitation water; RW: river water; SGW: shallow phreatic groundwater; DGW: deep confined groundwater; FLW: relatively fresh lake water; SLW: salt lake water. ^aOnly one representative sample data. ^bShallow phreatic water sample data.

Groundwater shows a similar hydrochemical evolution along the flow path. The average TDS values vary from 618 to 32 029 mg L⁻¹ for SGW and from 547 to 1401 mg L⁻¹ for DGW from the upstream area (Zone 2) to the middle to lower stream area (Zone 4). DGW is much fresher when contrasted with the SGW at the same location (Fig. 2c). There is essentially no difference in TDS between SGW and DGW from the central depression (Zone 5), with the values ranging from 336 229 to 361 200 mg L⁻¹ for SGW and 370 940 mg L⁻¹ for representative DGW (Table 1). Groundwater in the alluvial fan area (Zone 2) is dominated by ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, Cl⁻ and Na⁺with the concentrations ranging from 89 to 309 mg L⁻¹ for ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, from 90 to 437 mg L⁻¹ for Cl⁻ and from 79 to 232 mg L⁻¹ for Na⁺, respectively. To the middle to lower stream area (Zone 4), the dominant ions vary to Cl⁻ and Na⁺ for both SGW and DGW. The mean concentration of Cl⁻ is 11 550 mg L⁻¹ for SGW and 263 mg L⁻¹ for DGW, and the average concentration of Na⁺ is 10 464 mg L⁻¹ for SGW and 407 mg L⁻¹ for DGW. All groundwaters including SGW and DGW in the basin center (Zone 5) are dominated by Cl⁻, Na⁺ and Mg²⁺. SGW has the concentration ranging from 215 561 to 227 451 mg L⁻¹ for Cl⁻, from 12 388 to 35 713 mg L⁻¹ for Na⁺ and from 53 480 to 64 860 mg L⁻¹ for Mg²⁺. The concentration of representative DGW is 222 404 mg L⁻¹ for Cl⁻, 32 378 mg L⁻¹ for Na⁺ and 57 079 mg L⁻¹ for Mg²⁺. Overall, the water types of both SGW and DGW evolve from HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na type in the upstream area (Zone 2) to Cl-Na type in the middle to lower stream area (Zone 4), and eventually to Cl-Mg type in the low-lying depression (Zone 5) (Fig. 2b).

4.2 Stable and radio isotopes

The statistical summary of isotope results for precipitation, river water, lake water and groundwater can be found in Fig. 3 and Table 2. The representative snowmelt water in the Kunlun mountainous area (Zone 1) has a δD value of −77.0 ‰ and δ¹⁸O value of −11.9 ‰. The δD and δ¹⁸O values of precipitation in the mountainous area (Zone 1) are in the range of −85.3 ‰ to −71.6 ‰ and −10.9 ‰ to −9.3 ‰, with an average value of −75.2 ‰ and −10.0 ‰, respectively. The δD values of precipitation in the alluvial fan (Zone 2) range from −68.1 ‰ to −66.2 ‰, with the average value of −67.2 ‰ and δ¹⁸O values ranging from −10.1 ‰ to −9.7 ‰ with an average value of −9.9 ‰. The enrichment of stable hydrogen and oxygen isotopes in precipitation from the mountainous area to the basin reflects the secondary evaporation effect of precipitation in arid inland areas (Clark and Fritz, 1997). The δD and δ¹⁸O values of river water vary from −75.7 ‰ to −46.7 ‰ and between −11.1 ‰ and −4.8 ‰, respectively, showing a gradual enrichment trend along the river flow path. Fresh and salt lake waters are all significantly enriched in heavy isotopes with values ranging from −27.0 ‰ to −4.0 ‰ for δD and from 0.4 ‰ to 4.5 ‰ for δ¹⁸O. As shown in Fig. 3a and b, the SGW from the alluvial fan (Zone 2) to the middle to lower stream area (Zone 4) show an overall gradual enrichment trend along the flow path. In contrast, the DGW shows a significant depletion trend from the south to the north. Groundwater at different depths in the low-lying depression (Zone 5) are all brines with δD values ranging from −66.0 ‰ to −2.0 ‰ and the δ¹⁸O values ranging between −10.8 ‰ and −0.6 ‰, demonstrating relative enriched characteristics in contrast with the fresher groundwater in the upstream areas.

The ³H values range from 56.3 to 12.1 TU in the SGW and from 25.7 to <1 TU in the DGW along the groundwater flow path (Fig. 3d). The ¹⁴C activities in SGW vary from 57.9 to 11.9 pMC, and in DGW range from 49.15 to 0.7 pMC along the flow path (Fig. 3c). The spatial distributions of ³H and ¹⁴C results indicate increasing residence times for groundwaters in the aquifers from the south to north. While one shallow phreatic groundwater sample (G178) in the low-lying depression (Zone 5) was observed with relative high ³H content (18.9 TU), this may be caused by the rapid infiltration of surface water in the flood period.

Groundwater in the alluvial fan (Zone 2) has a high tritium content ranging from 20.0 to 56.3 TU with the average value of 35.5 TU, indicating recharge by modern water less than 60 years old. Shallow groundwater in the overflow zone (Zone 3) and middle to lower stream area (Zone 4) also has a relative high tritium content in the range of 12.1–25.7 TU with an average value of 21.1 TU for Zone 3 and 14.4–17.5 TU with an average value of 16.0 TU for Zone 4, and the representative shallow phreatic water adjacent to the salt lake (Zone 5) also shows a high tritium content of 18.9 TU, presenting modern water isotopic signatures. This may be caused by the mixture with the infiltrating modern surface water. DGW in the overflow zone (Zone 3) and middle to lower stream area (Zone 4) is with tritium content ranging from <1 to 10.1 TU for Zone 3 and <1 to 4.1 TU for Zone 4. The elevated tritium determined from several deep confined waters is most likely caused by mixtures with shallow phreatic water in the borehole; therefore, they cannot be used for groundwater age determination. Most DGW samples have a tritium content less than 1 TU, indicating they are not influenced by mixing with shallow groundwater in the boreholes. The age of this tritium-free DGW can be estimated using the radiocarbon activity.

Radiocarbon activity of groundwater can be significantly influenced by geochemical reactions (e.g., carbon minerals dissolution, isotopic exchange processes) during subsurface infiltration and in the aquifers (Cartwright et al., 2010b). It is therefore essential to correct the ¹⁴C activity on the total dissolved inorganic carbon (TDIC) before using it for groundwater age estimation. Many models, such as statistical models, geochemical models, and mixing models, were proposed for ¹⁴C activity correction. Most of the models are of limited interest due to the assumptions of fully closed systems or open systems, simplification, or even total ignorance of geochemical reactions beyond the recharge area. Using a model based on Carbon-13 is a good approach to correct the influence of geochemical reactions on ¹⁴C activity of TDIC, and suitable for both open and closed systems. The measured apparent ¹⁴C activity (¹⁴C_uncorr) of TDIC was corrected using δ¹³C as follows (Clark and Fritz, 1997):

$${}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{corr}}={}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{uncorr}}{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{rech}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{carb}}}{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{TDIC}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{carb}}}},$$

where ¹⁴C_corr is the corrected ¹⁴C activity of TDIC, δ¹³C_TDIC is the measured δ¹³C ratio of TDIC, δ¹³C_rech is the assumed initial δ¹³C ratio, and δ¹³C_carb is the δ¹³C ratio of carbonate being dissolved.

Groundwater in the study area is mainly recharged by Golmud River seepage in the upper alluvial fan located near parts of the Gobi desert where there is a lack of vegetation. The ¹⁴C activity and δ¹³C ratio of TDIC of the water would not be changed when infiltrating though the unsaturated zone. Thus, the δ¹³C_rech ratio should be equal or close to the atmospheric value (−6.4 ‰). δ¹³C_carb is close to 0 ‰ (Clark and Fritz, 1997). Only some of the tritium-free DGW samples in Zone 3 and 4 have measured δ¹³C data, and these were selected to calculate groundwater age using the aforementioned δ¹³C correction approach. The age of DGW in Zone 3 and  4 ranges from 2264 to 20 754 years along the flow paths. Due to the absence of radiocarbon data, the age of paleogroundwater in Zone 5 cannot be calculated, but it is certain that the age is more than 20 000 years, which was deduced from the oldest age of groundwater in Zone 4 (20 754 years).

4.3 Two-dimensional groundwater flow modeling

The groundwater flow model was calibrated using annual average hydraulic heads from 63 different shallow wells measured in 2015 along the cross section (not shown in figure). The calibration shows a good match between simulated and observed hydraulic heads, as demonstrated in Fig. 4. The comparison results show that the fit to observed hydraulic heads is better in the loess plain (including Zone 3, 4, 5) with the maximum deviation less than 0.8 m, while relative poor in the alluvial fan (Zone 2) with the maximum deviation less than 5 m. The deviation in the loess plain is mainly caused by the heterogeneity and anisotropy in lithology (Gu et al., 2017). The relative large deviation in results within the alluvial fan is most likely attributed to the steep hydraulic gradient (Islam et al., 2017) and larger seasonal fluctuation of hydraulic heads. Over the whole study area, the RMSE (root mean squared error) is only 1.57 m; therefore, the calibrated model can be used to reveal the groundwater flow pattern.

Figure 4 Comparison of observed and simulated hydraulic head values along the groundwater flow system, Golmud Watershed, China.

The estimated hydraulic parameters are shown in Table 3. The estimated values of K_h are 56.3 m d⁻¹ for gravel sand, 13.7 m d⁻¹ for sand, 0.62 m d⁻¹ for sandy silt, 0.13 m d⁻¹ for silt and 0.001 m d⁻¹ for clay. The anisotropy ratio of K_h∕K_v was estimated as 10 for gravel sand and sand, and 5 for sandy silt, silt and clay. These parameters are effective values under the assumption of homogeneity in each layer. The water budget analysis indicates a dynamic equilibrium state with the equilibrium difference of 0.62 %. Springs are the dominant discharge form, followed by evaporation and lake discharge, accounting for 76.81 %, 22.44 % and 1.37 %, respectively.

Table 3 Estimated parameters of different lithology from the Golmud Watershed, Qaidam Basin, China.

Figure 5 δ¹⁸O vs. δD diagram of precipitation, surface water and groundwater for the Golmud study area of the Qaidam Basin, China. (a) All data; (b) snow, precipitation and surface waters; (c) shallow phreatic waters; (d) deep confined waters.

5 Discussion

5.1 Water provenance and recharge characteristics

The δD and δ¹⁸O isotope analysis results for different water types are shown in Fig. 5a in relation to the global meteoric water line (GMWL: $\mathit{\delta }D=\mathrm{8}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{10}$) (Craig, 1961). The Golmud watershed local meteoric water line (LMWL: $\mathit{\delta }D=\mathrm{6.98}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{9.6}$) (Wang, 2014) and Golmud watershed local evaporation line (LEL: $\mathit{\delta }D=\mathrm{4.09}\times {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{28.1}$, R²=0.94), which is the linear regression line of river and lake water in the study area, are also shown in Fig. 5a. The slope and intercept of the LMWL (6.98 and 9.6) are lower than those of the GMWL (8 and 10) as a result of secondary evaporation that occurred during precipitation, reflecting the arid climatic characteristics of the study area (Dogramaci et al., 2012; Wang et al., 2017).

As shown in Fig. 5a, most of the surface water and groundwaters in the study area are situated close to the GMWL and LMWL, indicating a meteoric origin. However, the spatial distribution of precipitation is extremely uneven. Most of the precipitation occurred in the Kunlun mountainous area to the south. Precipitation in the basin is very limited (annual rainfall less than 50 mm) and in this area there is little effective recharge to the aquifers (Xiao et al., 2017). Thus, surface water and groundwater in the study area mainly originates from meteoric water (including precipitation and snowmelt) in the mountainous areas. River water (δD: −75.4 ‰ to −64.8 ‰, δ¹⁸O: −11.1 ‰ to −9.3 ‰) and groundwater (δD: −65.0 ‰, δ¹⁸O: −9.7 ‰) in the mountainous area (Zone 1) have similar stable water isotopic signatures to precipitation (δD: −85.3 ‰ to −71.6 ‰, δ¹⁸O: −10.9 ‰ to −9.3 ‰ ) and snowmelt water (δD: −77.0 ‰, δ¹⁸O: −11.9 ‰) values from the Kunlun mountainous area (Zone 1), indicating their direct recharge relationship (Fig. 5a). River waters flow towards the northern low-lying depression of the central basin and show a gradual enrichment trend due to intensive evaporation. Lake waters sampled from the low-lying depression (Zone 5) have the most enriched stable water isotope values and lie at the end of LEL defined by progressive evaporative enrichment of river water samples (Fig. 5b).

The δD and δ¹⁸O values of the SGW and DGW demonstrate different varying trends along the groundwater flow path. The SGW shows a gradual positive enrichment trend in heavy isotopes along the LEL (Fig. 5c), implying the influence of evaporation. For the alluvial fan (Zone 2), the δD and δ¹⁸O values of groundwater are very similar to that of river water in the alluvial fan (Zone 2) and groundwater in the mountainous area (Zone 1) (Table 2), indicating groundwater in the alluvial fan (Zone 2) is recharged directly by the seepage of river water and lateral inflow from the mountainous area, and out of the influence of evaporation. This is confirmed by similarities in major chemical composition (Fig. 2). The similar stable isotopic values also signify groundwater in the alluvial fan (Zone 2) has a short residence time, which is corroborated by elevated ³H (20.0–56.3 TU, mean value of 35.5 TU) indicating the residence time is less than 60 years based on ³H data (Xiao et al., 2017). SGW in the overflow zone (Zone 3) and the middle to lower stream area (Zone 4) has relative higher stable water isotope values compared with that in the alluvial fan and is plotted along the LEL, indicating SGW is influenced by evaporation from the overflow area (Zone 3) to the downstream. SGW in these two zones (Zone 3 and 4) also presents similar stable hydrogen and oxygen isotopic signatures as the river waters in the same area (Table 2), implying SGW has a very close hydraulic relationship with the rivers. The ³H values of the SGW in Zone 3 and 4 range from 12.1 to 25.7 TU and from 14.4 to 17.5 TU, respectively, with the mean value of 21.1 and 16.0 TU, suggesting that SGW in these two zones contains a large component of modern water or mixtures of old and modern water.

DGW in the overflow zone (Zone 3) and the middle to lower stream area (Zone 4) is observed to have a completely opposite evolution trend in that the δD and δ¹⁸O values become more depleted along the groundwater flow path (Fig. 5d). The depleted nature of the δD and δ¹⁸O values may have two interpretations: (1) these aquifers have another recharge region where rainfall with low δD and δ¹⁸O values occurs, or (2) the groundwater is ancient water recharged under colder climatic conditions (Chen et al., 2012; Awaleh et al., 2017). If (1) is the reason that groundwater would be more depleted in δD and δ¹⁸O along the groundwater flow paths, it is difficult to construct a mixing model that would supply more depleted waters along a deep flow path. According to the groundwater age estimated using ¹⁴C activity, the DGW in Zone 3 and 4 was recharged from 2264 BP to more than 20 754 BP (Holocene to late Pleistocene), which was a period when the climate changed from cold and wet condition (30 000 to 17 000 BP) to warm and dry conditions (14 000 BP to present) (Zhang et al., 2011). Consequently, it is believed that the depleted δD and δ¹⁸O waters in the DGW were recharged by paleowater under a colder climate relative to the present day. Similar findings were reported in the adjoining Nomhon watershed of the Qaidam Basin (Xiao et al., 2017). Additionally, this is consistent with the paleoclimate findings recorded using groundwater data from other basins of northwestern China (He et al., 2015; Huang et al., 2017).

Groundwaters in the low-lying depression area (Zone 5), regardless of depth, are all brines with TDS values greater than 100 000 mg L⁻¹. Given the tectonic activity and depocenter migration within the Qaidam Basin over geological history (Chen and Bowler, 1986; Zhang, 1987), groundwater in the low-lying depression area (Zone 5) has a large component of paleobrines migrated from western Qaidam Basin due to the uplift in the past (Huang and Han, 2007). According to the ¹⁴C age of DGW in Zone 4, the deduced age of DGW in the low-lying depression (Zone 5) is more than 20 000 years. SGW was observed with high ³H content (18.9 TU) adjacent to the terminal lake (G178) (Fig. 1 and Table 2), indicating mixing with leakage of modern surface water. As shown in Fig. 5a, most groundwaters in the basin center show a considerable deuterium excess, indicating that the original precipitation waters experienced considerable evaporation and vapor re-equilibration during recharge (Clark and Fritz, 1997).

5.2 Mechanisms controlling hydrochemistry

Figure 6 Diagrammatic representation showing the mechanisms controlling groundwater chemistry. (a) TDS vs. Na⁺ ∕ (Na⁺ + Ca²⁺); (b) TDS vs. Cl⁻ ∕ (Cl⁻ + ${\mathrm{HCO}}_{\mathrm{3}}^{-}$) (after Gibbs, 1970).

Figure 7 Bivariate plots (meq L⁻¹) of various ions in shallow phreatic and deep confined groundwater showed state (a) Na vs. Cl, (b) Na vs. (Cl + SO₄), (c) (Ca + Mg) vs. (HCO₃ + SO₄), (d) (Ca + Mg) vs. HCO₃, (e) Ca vs. SO₄, (f) (Ca + Mg)-(HCO₃ + SO₄) vs. [HCO₃ + SO₄].

Generally, the composition of natural groundwater is primarily controlled by the chemical composition of recharge waters, water–aquifer matrix interaction and groundwater residence time (Redwan and Moneim, 2015; Verma et al., 2016). As exhibited in the diagrams between TDS vs. Na⁺ ∕ (${\mathrm{Na}}^{+}+{\mathrm{Ca}}^{\mathrm{2}+}$) and Cl⁻ ∕ (${\mathrm{Cl}}^{-}+{\mathrm{HCO}}_{\mathrm{3}}^{-}$) (Fig. 6), the major mechanisms controlling groundwater chemistry are water–rock interaction and evaporation–mineral precipitation processes (Gibbs, 1970). Water–rock interaction processes dominate the controls on groundwater chemistry at all depths in the alluvial fan (Zone 2) due to the great depth and the negligible impact of evaporation. For the overflow zone (Zone 3) and the middle to lower stream area (Zone 4), the governing mechanisms for SGW change from water–rock interaction to evaporation–mineral precipitation due to the gradual decrease in groundwater depth and recharge inputs from waters having undergone the influence of intensive evaporation in that part of the basin. Nearly all DGW in this part of the flow system are controlled by water–rock interaction. Two DGW samples are observed to plot in the evaporation–crystallization domain (Fig. 6). This is due to a high TDS and oversaturation of evaporative minerals (such as aragonite, calcite and dolomite) in the groundwater resulting in mineral precipitation (crystallization). For the low-lying depression (Zone 5), evaporation has a significant influence on the chemistry of SGW, and crystallization (precipitation) of many mineral phases is the primary geochemical process controlling both the SGW and DGW chemistry.

In order to further constrain the sources of solute in groundwater, the relationships between various ions are compared (Fig. 7). The relation of Na⁺ vs. Cl⁻ shows both SGW and DGW from the piedmont to the middle to lower stream area (Zone 2, 3, 4) are plotted along the 1:1 line (Fig. 7a), suggesting that halite dissolution is potentially a primary process or source of Na⁺ and Cl⁻ mineralization in groundwater. The calculated results of halite saturation index (SI_halite<0) (Table 4) confirm that halite minerals of the aquifer matrix could be readily available to the groundwater. In addition, core drilling demonstrated that evaporated salts such as halite, calcium sulfate and sodium sulfate are widespread in the aquifer materials and can provide the solute source. Some of the SGW in Zone 4 is observed to have excess Na⁺ relative to Cl⁻ (Na ∕ Cl ratios equal to 1.2–3.8), while groundwater in Zone 5, regardless of the depth, shows deficiency of Na⁺ with respect to Cl⁻ (Na ∕ Cl ratios equal to 0.08–0.26), implying the existence of some other processes contributing Na⁺ not Cl⁻ to groundwater and changing the ratio of Na⁺ ∕ Cl⁻.

Table 4 Saturation index of selected minerals from the Golmud Watershed, Qaidam Basin, China.

One explanation for the excess of Na⁺ would be that the abundant Ca²⁺ and Mg²⁺ in fresh groundwater exchanges with the Na⁺ on the surface of clay minerals, which results in an increase in Na⁺ concentration and a decrease in Ca²⁺ and Mg²⁺ concentration in groundwater (Awaleh et al., 2017). The relationship of [$({\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{Mg}}^{\mathrm{2}+})$-$({\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-})$] vs. [$({\mathrm{Na}}^{+}+{\mathrm{K}}^{+})$-Cl⁻] (Fig. 7f) shows a regression line of $y=\mathrm{1.0016}x+\mathrm{4.9078}$ (R²=0.9966) and corroborates the contribution of cation exchange (Ca²⁺ or ${\mathrm{Mg}}^{\mathrm{2}+}+\mathrm{2}\mathrm{NaX}$ (solid) → $\mathrm{2}{\mathrm{Na}}^{+}+{\mathrm{CaX}}_{\mathrm{2}}$ or MgX₂ (solid)) (Verma et al., 2016). In addition, silicate weathering (e.g., 2NaAlSi₃O₈ (Albite) + 2CO₂+11H₂O → $\mathrm{2}{\mathrm{Na}}^{+}+{\mathrm{Al}}_{\mathrm{2}}{\mathrm{Si}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}(\mathrm{OH}{)}_{\mathrm{5}}$ (Kaolinite) + $\mathrm{3}{\mathrm{H}}_{\mathrm{4}}{\mathrm{SiO}}_{\mathrm{4}}+\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}$) in the aquifers could also contribute Na⁺, not Cl⁻, to groundwater (Guo et al., 2015). The ratio of Na⁺ ∕ (${\mathrm{Cl}}^{-}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) is around 1 (Fig. 7b), demonstrating mirabilite (Na₂SO₄⋅10H₂O) dissolution (${\mathrm{Na}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\cdot \mathrm{10}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to \mathrm{2}{\mathrm{Na}}^{+}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) is an additional strong possible process that could also be responsible for the excess of Na⁺ compared to Cl⁻ in groundwater (Jia et al., 2017). Groundwater in the basin with extremely high TDS concentration (more than 1 000 000 mg L⁻¹) has very low ratios of Na⁺ ∕ Cl⁻ (0.08–0.26) as a result of suspected reverse cation exchange (${\mathrm{Na}}^{+}+{\mathrm{CaX}}_{\mathrm{2}}$ or MgX₂ (solid) → $\mathrm{2}{\mathrm{Na}}^{+}+{\mathrm{CaX}}_{\mathrm{2}}$ or MgX₂ (solid)) (Fig. 7f).

The relationship between (${\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{Mg}}^{\mathrm{2}+}$) and (${\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) shows that almost all groundwater from the piedmont to the middle to lower stream area (Zone 2, 3, 4) are plotted along the 1:1 line (Fig. 7c), implying the dissolution of minerals such as gypsum, anhydrite, aragonite, calcite and dolomite are the potential ion sources to groundwater in the mineralization process (Dogramaci et al., 2012). As shown in Fig. 7d, nearly all groundwater data plot away from the equiline of (${\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{Mg}}^{\mathrm{2}+}$) vs. ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ (only three samples with the (${\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{Mg}}^{\mathrm{2}+}$) ∕ ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ ratio in the range of 0.8–1.2, 6 samples with the ratio range of 0.2–0.5, and others with the ratio range of 1.2–604.2), indicating that the Ca²⁺, Mg²⁺ and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ are not primarily derived from the dissolution of aragonite, calcite and dolomite. The saturation index values of aragonite, calcite and dolomite are all almost greater than 0 in all samples (Table 4), suggesting the dissolution of these three minerals must be minimal. The saturation index values of gypsum and anhydrite for groundwater in these areas are all below zero (Table 4), corroborating the contribution of gypsum and anhydrite dissolution for groundwater mineralization. The deficiency of Ca²⁺ compared to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ of groundwater (73.5 % of samples with the Ca ∕ SO₄ ratio less than 0.8) presented in Fig. 7e is most likely as a result of the aforementioned mirabilite (Na₂SO₄⋅10H₂O) dissolution and cation exchange. As mentioned earlier, the redox conditions of the deep confined aquifers in Zone 4 have evolved to a reduced environment, but due to the extremely low organic carbon content in the sediments (Bowler et al., 1986; Chen and Bowler, 1986), sulfate reduction has a very limited influence on groundwater chemical evolution. This is also the reason that groundwater in the downstream area (Zone 4 and 5) has an abundant content of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in contrast to Ca²⁺.

Groundwater in the low-lying depression (Zone 5) has extremely high TDS values (>300 000 mg L⁻¹) (Table 1) and almost all minerals are oversaturated (SI > 0) (Table 4); therefore, precipitation (crystallization) of minerals is the primary geochemical process in this part of the aquifer (Li et al., 2010). In addition, reverse cation exchange interaction and evaporation, which can be confirmed by the relationship of [$({\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{Mg}}^{\mathrm{2}+})$-$({\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-})$] vs. [$({\mathrm{Na}}^{+}+{\mathrm{K}}^{+})$-Cl⁻] (Fig. 7f) and the relation of stable water isotopes (Fig. 4), respectively, are also important mechanisms governing groundwater chemistry. Surface water has significant influences on the geochemical processes that occurred in the shallow aquifers. In the wet season, a large amount of relatively freshwater can reach the low-lying depression area (Zone 5) and infiltrate to the shallow aquifers. This would dilute the groundwater and dissolve the evaporated salts in the aquifers.

Figure 8 Conceptual model of groundwater flow and hydrochemical evolution in the Golmud watershed, China.

5.3 Groundwater flow and hydrogeochemical evolution

Theoretically, three types of groundwater flow systems, namely local, intermediate and regional, may occur in a large basin, and each flow system has its own characteristics based on aspects of flow path, recharge origin, cycle depth, cycle amount, residence time, discharge position, hydrochemistry and controlling mechanisms (Toth, 1963). The cross-sectional groundwater flow modeling results demonstrated the groundwater flow paths in the study area are strictly controlled by distribution of the lithology (Fig. 8). Groundwater flow lines are shown to be upward convex in shape at the front of the alluvial fan and the middle-stream area due to an increase in relatively poor permeability due to the addition of finer, less permeable stratigraphic material. Based on the distribution of flow lines, three groundwater flow systems including local, intermediate and regional system were identified in the study area (Fig. 8).

The local groundwater flow system occurs in the shallow part of the alluvial fan (Zone 2) and overflow zone (Zone 3) with the deepest cycle depth within 250 m of the surface. This flow system obtains recharge water along the Golmud River flow path and discharges at the overflow zone. The water cycle quantity of the local system estimated by modeling accounts for approximately 83 % of the total quantity of groundwater in the watershed. Groundwater chemistry is mainly controlled by water–rock interaction and there appears to be very little evaporation. Groundwater has a rapid velocity in this part of the system, with a residence time of less than 60 years. As a result, groundwater here is fresh, with TDS values less than 1000 mg L⁻¹, and the water type is mainly HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na. This system is the main source of water supply for Golmud city.

The intermediate flow system occurs below the local system and is recharged by river water seepage near the upper part of the alluvial fan. Groundwater flows to lower elevations towards the north and reaches its deepest cycle depth near 600 m at the middle part of the alluvial fan (Zone 2). Due to the increase in aquitards, water flow lines are presented as upward convex shapes at the middle to lower part of alluvial fan (Zone 2). Groundwater flow is constrained by two continuous aquitards (clay layers) at depths of 60 and 290 m (Fig. 8), respectively, at the front of the alluvial fan and overflow zone. The intermediate flow system discharges between the lower overflow zone (Zone 3) and the middle to lower stream area (Zone 4), as evidenced by springs and surface evaporation. The total cycle water quantity of the intermediate system accounts for approximately 14 % of the total cycle groundwater amount in the watershed. Aquifers of this system in the alluvial fan have higher renewal rates due to their increased permeability, compared to those in the lower overflow zone (Zone 3) and the middle to lower stream area (Zone 4) that have relative low renewal rates as a result of a lithology dominated by finer sediments, with groundwater residence times of about 4000 years. Hydrochemistry is dominantly controlled by water–rock interaction, and also strongly influenced by evaporation within the discharge area. Because of the short residence time and a shortage of chemical solutes in the aquifer material, groundwater in the alluvial fan (Zone 2) generally maintains recharge water chemical characteristics that are fresh and of HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na type. Sufficient solutes in the aquifer medium and intensive evaporation in the fine soil plain results in the groundwaters gradually evolving to be brackish water and in some cases saline waters.

The regional groundwater flow system occurs under the intermediate system and is recharged at the upper part of the alluvial fan by river water seepage and lateral flow within the mountainous area, and it discharges at the basin center into terminal lakes resulting in evaporation. Groundwater flow paths are significantly controlled by the lithology (Fig. 8), and divided from the intermediate system by a continuous aquitard at a depth of 290 m. Aquifers of this system have very low water renewal rates with residence times up to and greater than 20 000 years. The modeled cycle water quantity of the regional system is only approximately 3 %. Groundwater chemistry is mainly influenced by water–rock interaction in this system, except for shallow aquifers in the discharge area (Zone 5), which are strongly influenced by evaporation. Due to the substantial difference in residence time, water–rock interaction results in much different hydrochemical characteristics from the other aquifer systems. Groundwater that was presented as freshwater with a dominant water type of HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na in the alluvial fan (Zone 2) and overflow zone (Zone 3) become a brackish water type (HCO₃ ⋅ Cl-Na) and a saline water type (Cl-Na) in the middle to lower stream area (Zone 4) and has evolved to be a brine water type mainly composed of Cl-Mg in the low-lying discharge area (Zone 5).

6 Conclusions

Previous studies on arid closed basins such as the Great Artesian Basin, Murray Basin, Death Valley and Minqin Basin have established a lot of typical groundwater circulation and evolution regimes. However, the Qaidam basin, a typical arid sedimentary closed basin formed with the uplift of the Tibetan plateau, has groundwater circulation patterns characterized by the complex tectonic activities, paleoclimate variation, arid climate characteristics, sedimentary lithology, and systematic evolution from fresh- to saltwater. Studies of this basin can enhance the understanding of groundwater origin, flow regime and hydrogeochemical evolution in such complex tectonic influenced arid sedimentary closed basins worldwide. Integration of hydrogeochemistry, isotopes and two-dimensional groundwater flow modeling was used to obtain insight into the hydrogeology in a typical arid endorheic watershed represented by the Qaidam Basin, Tibetan plateau. A number of key findings have come out of this study.

The groundwater in the basin originates from precipitation and meltwater in the mountainous areas to the south. Groundwater in the alluvial fan is recharged directly as a result of modern river water seepage and mountainous lateral inflow and has a rapid flow rate. Shallow phreatic waters in the overflow zone and the middle to lower stream area are supported by local and intermediate groundwater flow systems and have a close chemical and isotopic relationship with surface water. Deep confined groundwater in the overflow zone and the middle to lower stream area is recharged from paleometeoric water during the late Pleistocene and Holocene under a cold climate based on the results for stable water isotopic analyses. Groundwater in the low-lying depression (basin center) is made up of ancient brines that have possibly migrated from the western Qaidam Basin due to the uplift of the western basin in the geological past. Shallow phreatic aquifers in the low-lying depression (basin center) are also seasonally recharged by modern surface water during flooding periods.

Groundwater in the study area evolves from freshwater to brine water along the flow path. The hydrochemistry of groundwater in the alluvial fan is dominantly controlled by mineral dissolution and cation exchange, and occurs as slightly alkaline water with TDS values less than 1000 mg L⁻¹ and a water type with a composition of HCO₃ ⋅ Cl-Ca ⋅ Mg ⋅ Na. Deep confined groundwater chemistry in the overflow zone and middle to lower stream area is also controlled by mineral dissolution and cation exchange, as a result of longer residence times in the aquifers and shows a trend evolving from freshwater to brackish water and finally saline water with increasing solute inputs along the flow paths. As well as water–rock interaction, shallow phreatic water is also affected by intensive evaporation, and therefore, these waters can be much saltier than deep confined water. Groundwater in the low-lying depression (basin center) is composed of brine water, and the mineral precipitation coupled with reverse cation exchange is the dominant geochemical process controlling water chemistry. The impact of evaporation is also one of the important geochemical processes in the shallow phreatic aquifers, which can accelerate evaporated mineral precipitation. Salt dissolution occasionally occurred in the low-lying depression (basin center) during flood periods due to the infiltration of large amounts of fresh surface water.

Three different hierarchical groundwater flow systems were identified using the cross-sectional model. The continuous aquitards at depths of 60, 290 and 450 m have significant constraints on groundwater flow. The local flow system occurs in the shallow part of the alluvial fan and overflow zone and discharges in the overflow zone, with the deepest cycle depth within 250 m of surface. The intermediate system occurs below the local system and discharges in the lower overflow zone and middle to lower stream area with the deepest cycle depth reaching 600 m below surface. The regional system was separated from the intermediate system by a continuous aquitard at a depth of 290 m and discharges in the low-lying depression (basin center). Our calculation shows that the discharge water quantity of these three systems accounts for approximately 83 %, 14 % and 3 %, respectively.

This study enhanced the understanding of the origin, flow pattern, hydrochemical evolution and controlling mechanisms of the regional groundwater systems in the Qaidam Basin. These results can provide fundamental information for coping with future issues such as water conflicts, salt lake exploitation and climate warming in the basin and also provide references for understanding the hydrogeological processes in other similar endorheic watersheds of northwestern China and elsewhere in the world.