Saline groundwater evolution in the Luanhe River delta (China) during the Holocene: hydrochemical, isotopic, and sedimentary evidence

Since the Quaternary Period, paleo-seawater intrusions have been suggested to explain the observed saline groundwater that extends far inland in coastal zones. The Luanhe River delta (northwest coast of the Bohai Sea, China) is characterized by the distribution of saline, brine, brackish, and fresh groundwater from the coastline inland. The groundwater in this region exhibits a wide range of total dissolved solids (TDS): 0.38–125.9 g L⁻¹. Meanwhile, previous studies have revealed that this area was significantly affected by Holocene marine transgression. This study used hydrochemical, isotopic, and sedimentological methods to investigate groundwater salinization processes in the Luanhe River delta and its links to paleo-environmental settings. Isotopic results (²H, ¹⁸O, ¹⁴C) allowed old groundwater recharge to be distinguished from new groundwater recharge. Hydrochemical analysis using the PHREEQC code indicated that the salt in saline and brine groundwater originates from a marine source. The ¹⁸O–Cl relationship diagram yields three-end-member groundwater mixing, and two mixing scenarios are suggested to explain the freshening and salinization processes in the study area. When this was interpreted along with data from paleo-environmental sediments, we found that groundwater salinization may have occurred since the Holocene marine transgression. The brine is characterized by radiocarbon activities of ∼ 50–85 pMC and relatively depleted stable isotopes, which are associated with seawater evaporation in the ancient lagoon during delta progradation and mixing with deeper fresh groundwater, which was probably recharged in the cold Late Pleistocene. The brackish and fresh groundwaters are characterized by river-like stable isotope values, where high radiocarbon activities (74.3–105.9 pMC) were formed after the washing out of the salinized aquifer by surface water in the delta plain. This study presents an approach that utilizes geochemical indicator analysis with paleo-geographic reconstruction to better assess groundwater evolutionary patterns in coastal aquifers.

1 Introduction

It is estimated that around 40 % of the world's population live in coastal areas (UN Atlas, 2010). Groundwater is an important freshwater resource for domestic consumption and agricultural activities in this region (Cary et al., 2015; Jayathunga et al., 2020). However, groundwater salinization poses a significant threat to everyday living and development activities (Tulipano, 2005; de Montety et al., 2008). In recent decades, groundwater salinization in coastal zones has prompted widespread concern and has been widely studied. On the one hand, seawater intrusion due to groundwater pumping is an important salinization process in coastal aquifers (Reilly and Goodman, 1985; Werner, 2010, 2013; Han and Currell, 2018). On the other hand, groundwater salinization caused by paleo-seawater intrusion in response to Quaternary changes in the global sea level has been reported in many coastal zones worldwide (Edmunds, 2001; Akouvi, 2008; Santucci et al., 2016, Larsen et al., 2017).

Coastal aquifers are linked to the ocean and continental hydrological cycles (Ferguson and Gleeson, 2012), both of which are influenced by natural and human-induced change (Jiao and Post, 2019). In the natural state, there is a steady-state seawater–freshwater interface that extends inland from the coastal line (Costall et al., 2020). During the Quaternary Period, however, sea-level fluctuations on geological timescales have caused the interface to change, allowing seawater intrusion during transgression events and freshwater flushing during glacial low-sea-level periods, as is evident from the hydrochemical characteristics of groundwater in coastal aquifers (Kooi et al., 2000; Sanford, 2010; Aquilina et al., 2015; Lee et al., 2016). In addition, the hypersaline groundwater found in coastal zones, particularly brine groundwater with a salinity 2–4 times that of seawater, cannot be explained solely using a seawater intrusion model (Sola et al., 2014; Han et al., 2020), and paleo-environment settings must be taken into consideration (van Engelen et al., 2019). Some studies, for example, attribute the presence of brine in Mediterranean countries to the evaporation of seawater in the lagoon system during the Holocene transgression (Giambastiani et al., 2013; Vallejos et al., 2018).

The Bohai Sea of northern China was affected by Late Pleistocene transgressive–regressive cycles, which caused paleo-saltwater intrusions with various salinities along the coastal aquifers (Du et al., 2015; Li et al., 2017). Several studies have applied geochemical methods to elucidate the origin of saline groundwater and salinization processes with an anthropogenic influence, including the induced mixing of brine water from adjacent aquifers caused by groundwater overexploitation in Laizhou Bay (Han et al., 2011, 2014; Liu et al., 2017; Qi et al., 2019). However, the association between groundwater salinization (especially brine formation) and paleo-environmental implications are still not clear. Thus, this study applied a range of chemical, isotopic, and sedimentary indicators to examine the groundwater salinization processes in the Luanhe River delta (situated along the northwestern coast of the Bohai Sea) in relation to recharge, salt sources, mixing behavior, and paleo-geographic evolution. The overall goal was to understand the evolutionary pattern of groundwater as influenced by transgression/regression events on a geologic timescale. Our findings are relevant to aquifer remediation activities in the region as well as other similar sedimentary environments around the world.

2 Background on the study area

The study area is located in northeastern Hebei Province, China, on the west coast of Bohai (Fig. 1a). It consists of an alluvial fan and coastal delta bounded by a Holocene maximum transgression line (Xue, 2016). The delta area can be further divided into two parts: the old delta between the Douhe River and the Suhe River and the new delta between the Suhe River and the modern Luanhe River (He et al., 2020). In terms of geomorphology, the study area is inclined to the south and southwest with a slope of about 0.04 ‰–2 ‰. The temperate monsoon climate of the area affects its average annual temperature of 12.5 ^∘C and annual rainfall of 601 mm (1956–2010), with 80 % of the annual rainfall occurring between July and September.

2.1 Hydrogeology

The thickness of Quaternary sediments in the study area is about 400–500 m. According to their lithology and hydrogeological characteristics, there are four distinct Quaternary aquifers (Fig. 1b). The first Holocene aquifer (Q₄) is a phreatic or semiconfined aquifer with a bottom depth of 15–30 m, is primarily composed of fine sand and slit, and involves fresh, brackish, saline, and brine groundwater (Dang et al., 2020). The second Late Pleistocene aquifer (Q₃), the third Middle Pleistocene aquifer (Q₂), and the fourth Early Pleistocene aquifer (Q₁) have bottom depths of 120–170, 250–350, and 350–550 m, respectively. They are confined aquifers primarily made of medium sand and gravel (Niu et al., 2019). The first aquifer is mainly recharged by meteoric precipitation and lateral infiltration of surface water (Li et al., 2013). The groundwater from the first aquifer is widely extracted for irrigation in the alluvial fan area. The largest salt farm in North China, the Daqinghe Salt Farm, uses shallow brine groundwater for salt production in the delta area, where there is only minor agricultural activity. Except in the alluvial fan area, the circulation between phreatic and confined aquifers is weak. The deep groundwater in the second, third, and fourth aquifers is mainly recharged by a surrounding mountain range and mainly discharged by human pumping (Ma et al., 2014).

Figure 1 (a) Location map of the study area. Also shown are the locations of the sampling site and published cores in the Luanhe River delta. Cores LT05, HCZ, BXZK01, BXZK02 and BXZK03 are from He et al. (2020); cores NP05, NP03, DY01, DQ03, DQ04, DY02, MT04, BG07, FB01, A02, and TH04 are from Xu et al. (2020); core LQZ04 is from Cheng et al. (2020); core FG01 is from Xu et al. (2011); core Bai03 is from Li and Wang (1983); and core HCZ is from Peng et al. (1981). (b) Hydrogeological cross section (A–A' in Fig. 1a) of the study area; modified from Ma et al. (2014).

2.2 Sedimentary evolution since the Late Pleistocene

Previous studies have shown that in the study region, the salt–fresh groundwater interface gradually deepens from land (depth of ∼ 5 m) to sea (depth of ∼ 100 m), as shown in Fig. 1b, with salt groundwater primarily occurring in the first aquifer of the delta area (Li et al., 2013; Ma et al., 2014). According to a stratigraphic transect along the present coastline (Fig. 2), the series stratigraphic architecture of the first aquifer consists of Late Pleistocene continental facies – Holocene marine facies – Holocene delta facies – modern continental facies or artificial fill, indicating that the sediments of the first aquifer have been deposited through lowstand continental accumulation, marine transgression, and highstand progradation since the Late Pleistocene.

The seawater did not reach the modern coastline from the Last Glacial Maximum to the early Holocene (about 30–9 ka BP). The Luanhe alluvial fan was active during this period (He et al., 2020). Since about 9000 a BP, the Holocene marine transgression has approached the present coastline (Xu et al., 2020), and Holocene marine sediments developed due to the sea-level rise from 9–7 ka BP. The Holocene marine transgression had reached its maximum inland area 20 km from the modern coastline by about 7 ka BP (Gao et al., 1981; Peng et al., 1981; Xue, 2014, 2016) (see the Holocene transgression line in Fig. 1). The modern coastal plain was formed through the accumulation of the highstand prograding delta on top of Holocene marine strata, together with artificial fill. In addition, lagoons are important components of the Luanhe River delta (Feng and Zhang, 1998). Based on the records of lagoon facies in the published cores for this region, the approximate distribution range of the buried lagoon is shown as a purple dashed line in Fig. 1a.

Figure 2 Stratigraphic transect along the present coastline of the Luanhe River delta; modified from He et al. (2020).

3 Methods

In total, 45 water samples were collected from the Luanhe River delta, including 38 groundwater samples, five surface water samples, one local rain water sample, and one Bohai seawater sample, during four sampling campaigns from October 2016 to June 2020. Groundwater samples were divided into shallow groundwater samples and deep groundwater samples, which were pumped from unconfined and confined aquifers, respectively. The surface water included two Suhe River water samples and two Luanhe River water samples. Due to the artificial fill that has modified the coastal landscape, it was difficult to locate the modern lagoon environment. However, during the investigation, it was found that the Daqinghe Salt Farm in this area extracts seawater into the evaporation pond. The mixture of seawater and meteoric water is subject to evaporation, leading to the formation of concentrated saline water (CSW) in the pond, similar to the formation of CSW in a coastal lagoon (Stumpp et al., 2014). Thus, one CSW sample (sample P18) from the evaporation pond was collected.

Water types were classified according to Zhou (2013) into freshwater (TDS < 1 g L⁻¹), brackish water (TDS = 1–3 g L⁻¹), saline water (TDS = 3–50 g L⁻¹), and brine (TDS > 50 g L⁻¹). Groundwater sampling depths and pH values were measured on site using Eureka Manta+(water-quality multiprobe sonde). The concentrations of K⁺, Na⁺, Ca²⁺, Mg²⁺, and Br⁻ ions were measured using inductively coupled plasma analysis (ICAP-7400), while SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and Cl⁻ ions were determined using ion chromatography (ICS-600). The HCO${}_{\mathrm{3}}^{-}$ concentrations of samples were measured using titration. The resulting hydrochemical data are listed in Table S1 (see the Supplement). The stable isotope concentrations (δD, δ¹⁸O) in the water samples G02-10, G06-10, G03-05, G04-40, G05-10, G05-46, G07-27, P07-20, P08-30, P09-30, P10-30, P11-20, P12-40 P14-15, P07-100, P13-200, P14-300, P15-150, P16-100, P17-200, P18, LH01, LH02, SH01, SH02, SW01, and R1 were tested at the Experimental & Testing Center of Marine Geology, Ministry of Natural Resources, China, using high-temperature pyrolysis–isotope ratio mass spectrometry. The values of δ¹⁸O and δD were calculated with respect to Vienna Standard Mean Ocean Water (VSMOW), and the uncertainties in δD and δ¹⁸O were ±1.0 ‰ and ±0.2 ‰, respectively. The radioisotope measurements (AMS ¹⁴C) of the groundwater samples P14-300, P15-150, and P16-100 were performed at the Pilot National Laboratory for Marine Science and Technology. Analysis of the stable isotopes (δD, δ¹⁸O, ¹³C) and the radioisotope in groundwater samples G10-10, G03-20, G04-15, G05-30, G06-15, G07-15, G08-15, G08-40, G09-15, G09-40, G10-10, and G10-30 was performed at the Beta Analytic testing laboratory, where the δ¹⁸O and δD values were also calculated with respect to VSMOW, and the uncertainties in δD and δ¹⁸O are listed in Table S1. The ¹⁴C age of groundwater was calculated using the following equation (Clark and Fritz, 1997): $bt=-\mathrm{8267}\cdot \ln \left(a_{\mathrm{t}}^{\mathrm{14}}\mathrm{C}/q\cdot a_{\mathrm{0}}^{\mathrm{14}}\mathrm{C}\right)$, where t is the radiocarbon age in years Before Present (a BP); $a_t^{\mathrm{14}}$C is the measured ¹⁴C activity in % of modern carbon (pMC); $a_{\mathrm{0}}^{\mathrm{14}}$C is the derived modern ¹⁴C activity of soil; and q is a corrective factor accounting for the dissolution of calcite, which is assumed to be free of ¹⁴C and therefore dilutes the initial ¹⁴C activity of aqueous dissolved inorganic carbon (DIC) in recharged water. The results for ¹³C and ¹⁴C and the uncorrected residence times are listed in Table S2.

4 Results

4.1 Hydrochemistry

Except for P13-200 (TDS = 1.617 g L⁻¹, indicating brackish water), all the deep groundwater samples in the study region are freshwater. Deep groundwater hydrochemical forms shift from Ca-HCO₃ to Na-HCO₃ upon moving from the land to the sea (Fig. 3). For the shallow aquifer, the horizontal salt–fresh groundwater interface corresponds more closely to the maximum Holocene transgression line (see Fig. 1a). The Ca-HCO₃ type of shallow fresh groundwater is primarily distributed in the alluvial fan region. The brackish and low-TDS saline groundwaters, which vary from Ca-HCO₃ to Na-HCO₃ and Na-Cl types, occur mainly in the upper aquifer (depths of 0–15 m) of the delta area, while the lower part (depths of 20–40 m) contains Na-Cl-type saline and brine groundwaters with high TDS. Moreover, in terms of the horizontal distribution of salinity, the groundwater TDS tends to decrease from west to east, with the TDS values of saline and brine groundwater generally ranging from 16.57 to 125.97 g L⁻¹ in the old delta (western delta) and from 3.26 to 52.48 g L⁻¹ in the new delta (eastern delta).

Figure 3 Piper diagram of the various water samples.

4.2 The stable isotopes ²H and ¹⁸O

Figure 4 shows the relationship between deuterium and oxygen-18. The global meteoric water line (GMWL, δ²H = 8 ⋅ δ¹⁸O + 10) is from Craig (1961), while the local meteoric water line (LMWL, δ²H = 6.6 ⋅ δ¹⁸O + 0.3) is based on δ²H and δ¹⁸O isotope data (1985–2003, mean monthly rainfall values) from the Tianjin station about 100 km southwest of the study area (IAEA/WMO, 2006). The deep groundwater samples exhibit depleted values of stable isotopes, with values of δ²H ranging from −75.52 ‰ to −57.06 ‰ and those of δ¹⁸O from −9.82 ‰ to −7.61 ‰. Shallow groundwater samples have higher hydrogen and oxygen isotope levels, ranging from −64.6 % to −22.46 % for δ²H and −8.74 % to −2.07 % for δ¹⁸O. While the relatively small overall values of fresh and brackish groundwater samples are similar to those of the river samples, saline and brine groundwater samples generally plotted below the LMWL or GMWL, meaning that the water was subjected to evaporation prior to becoming groundwater recharge (Gibson et al., 1993), or that multiple-end-member mixing processes were involved (Han et al., 2011).

Figure 4 Stable isotope compositions of different water samples. Seawater mixing line: mixing between deep fresh groundwater and seawater; CSW mixing line: mixing between deep fresh groundwater and CSW.

4.3 Groundwater residence times

The measured ¹⁴C activities of groundwater samples range from 0.774 to 105.9 pMC (Table S2). A plot of sampling depth versus ¹⁴C is shown in Fig. 5; this elucidates a negative correlation, showing that the variation in ¹⁴C activity can be attributed to radioactive decay in the aquifer. There are multiple processes that can impact the ¹⁴C properties, including groundwater mixing and dispersion, long-term variation in atmospheric ¹⁴C, and free ¹⁴C dilution (e.g., carbonate dissolution) (Cartwright et al., 2020). Because the relative impacts of these processes are not well established in the study area, the uncertainty regarding the correction of radiocarbon ages to real groundwater ages is very high. Consequently, we estimate the groundwater age as a residence time range. The uncorrected age is considered the maximum age, while the corrected age is considered the minimum age. Corrected ages are determined based on two hypothetical models of carbonate dissolution, which is the main influence on the ¹⁴C contents of water samples (Lee et al., 2016).

Figure 5 shows that the ¹⁴C activities in the shallow groundwater samples are within 30.6–105.9 pMC. These values indicate a relatively modern recharge before the atmospheric nuclear testing period of the 1950s and 1960s. The radiocarbon activities in the deep fresh groundwater samples are less than 12 pMC, which is consistent with the paleo-water recharge. This indicates that there are weak connections between the shallow and deep aquifers. Therefore, we assume that the shallow aquifer is an open system and the deep aquifer is a closed system. The δ¹³C mixing and chemical mass balance (CMB) models are used to estimate the corrective factor q (Clark and Fritz, 1997).

For the δ¹³C mixing model (Pearson and Hanshaw, 1970), $q=(\mathit{\delta }{}^{\mathrm{13}}{C}_{\mathrm{DIC}}-\mathit{\delta }{}^{\mathrm{13}}{C}_{\mathrm{CARB}})/(\mathit{\delta }{}^{\mathrm{13}}{C}_{\mathrm{RECH}}-\mathit{\delta }{}^{\mathrm{13}}{C}_{\mathrm{CARB}})$, where δ¹³C_DIC is the measured δ¹³C of DIC in groundwater; δ¹³C_CARB is the δ¹³C of DIC from dissolved soil minerals (δ¹³C_CARB= 1.5 ‰ from Chen et al., 2003); and δ¹³C_RECH is the δ¹³C in water when it reaches the saturation zone. In this study, we use a δ¹³C_RECH of −15 ‰, which has been suggested as appropriate for soils dominated by C₄ plants in northern China (Currell et al., 2010). This model yields some relatively low q values (0.59 for G06-15 and 0.65 for G08-15), possibly because several unaccounted-for factors contribute to variable δ¹³C_RECH values, e.g., local methanogenesis and pH or temperatures in the soil zones.

For the CMB model, $q={\mathrm{mDIC}}_{\mathrm{rech}}/{\mathrm{mDIC}}_{\mathrm{final}}$, where mDIC_rech is the DIC molar concentration in the recharge water and mDIC_final is the DIC molar concentration in the final groundwater. mDIC_final was calculated using ${\mathrm{mDIC}}_{\mathrm{final}}={\mathrm{mDIC}}_{\mathrm{rech}}+[\mathrm{mCa}+{\text{Mg-SO}}_{\mathrm{4}}+\mathrm{0.5}(\mathrm{Na}+\text{K-Cl}\left)\right]$ (Fontes and Garnier, 1979). DIC_rech was mainly HCO₃ in the recharge water when pH values were between 6.4 and 10.3, and the carbonate equilibrium constant varied with temperature (Clark and Fritz, 1997). mDIC_rech was calculated from the estimated pH and temperature conditions in the recharge environment; e.g., at pH = 6 and T = 15 ^∘C, mDIC_rech= 10 mmol L⁻¹ (Currell et al., 2010).

The corrected radiocarbon ages are shown in Table S2. The residence time of deep groundwater ranges from 15 959 to 39 050 a BP, which is significantly longer than that of the groundwater in the shallow aquifer (9510 a BP to modern). Moreover, most brackish and fresh groundwater ages are modern, while brine has a longer residence period (5590–1245 a BP) and a broad range of residence times.

Figure 5 ¹⁴C activity with sampling depth in groundwater.

5 Discussion

5.1 Isotopic analysis of groundwater origin and recharge

Deuterium and oxygen-18 are good tracers of groundwater origin and climatic conditions during recharge periods (Clark and Fritz, 1997). When combined with groundwater residence times, they can also identify modern recharge and paleo-recharge (Han et al., 2014).

The depletion of ¹⁸O and ²H values in the deep fresh groundwater (Fig. 4) can be attributed to a cold climate (Kreuzer et al., 2009) and the residence times (33 951–39 050 a BP) of the P15-150 and P14-300 samples, which may suggest that there was a recharge during the last glacial maximum. The stable isotopes in P16-100 are heavier, reflecting a recharge history linked to the warm climate in the previous deglaciation (Hendry and Wassenaar, 2000). The stable isotope values of river samples are similar to those of the shallow brackish and fresh groundwater samples, which are approximately modern, indicating lateral recharge of surface water locally. Meanwhile, in Fig. 4, G03-5 is close to the rainfall sample, indicating that modern precipitation is a new recharge source. The trend towards δ²H and δ¹⁸O enrichment in brine and saline groundwater can be attributed to seawater infiltration during the Holocene transgression period, which has been confirmed by other studies of the Bohai Sea coast (Li et al., 2017; Du et al., 2016). Additionally, due to the mixing of meteoric water and the subsequent nonequilibrium fractionation of hydrogen isotopes during evaporation (Clark and Fritz, 1997), the CSW sample is characterized by ¹⁸O enrichment but ²H depletion compared to seawater.

5.2 Hydrochemical analysis of sources of salinity

To distinguish the sources of groundwater salinity, the PHREEQC code (Parkhurst and Appelo, 2013) was used to measure and plot the theoretical seawater–freshwater mixing line (“mixing line”) and seawater evaporation line (“evaporation line”) using hydrogeochemical modeling. The two simulation effects serve as references for groundwater hydrochemical characteristics (Figs. 6 and 7). For the Na-Cl (Fig. 6a), Mg-Cl (Fig. 6b), and Br-Cl (Fig. 7a) diagrams, the measured brackish, saline, and brine groundwater samples fit quite well to modeling mixing lines, and evaporation lines follow a linear trend from the least to the most saline. This strongly demonstrates that the salt in these water samples is mainly of marine origin. The major ion concentrations in some samples (such as brine) are higher than those in seawater, suggesting that the enriched ions are associated with evaporation processes rather than seawater intrusion (Colombani et al, 2017).

Moreover, the samples deviate from the modeling lines (Fig. 6c and d), indicating that there may be other hydrogeochemical processes that are responsible for the modified ionic compositions (Giambastiani et al., 2013). (1) Ca²⁺ depletion is seen for the P18 and P12 samples in Fig. 6d. This phenomenon is likely explained by gypsum (CaSO₄) precipitation. The evaporation line reveals that the Ca²⁺ composition of evaporating seawater follows a hooked trajectory (Fig. 6d). During evaporation to the point of gypsum saturation, the Ca²⁺ concentration of the residual CSW progressively decreases. (2) Ca²⁺ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ are present in excess in most of the fresh and brackish samples (Fig. 6c and d), which can be attributed to mineral dissolution along with stream water recharging (such as gypsum dissolution), highlighting the occurrence of some degree of dilution with continental runoff since the Holocene regression. (3) Decomposition of the organic matter abundant in marine or lagoon facies sediments can result in the release of bromide ions, thus pushing the Br $/$ Cl ratios of saline groundwater samples higher than the mixing line (Fig. 7b).

Figure 6 Hydrochemical relationships between Cl and major ions of the measured samples and simulated results in groundwater. Seawater–freshwater mixing line: theoretical mixing between seawater and deep fresh groundwater. Blue numbers: mixing ratios of seawater. Seawater evaporation line: theoretical evaporation of Bohai seawater. Red numbers: evaporation rates. G: precipitation of gypsum. H: precipitation of halite. Dissolution: possible gypsum dissolution along with river water recharging. The symbols used for the samples are the same as in Fig. 6.

Figure 7 Relationship between the chloride and bromide contents in water samples. Symbols used are the same as in Fig. 6.

5.3 Mixing processes

Figure 8 depicts the relationship between δ¹⁸O and Cl⁻ in different water samples. There are higher Cl⁻ concentrations and lower δ¹⁸O values in brine samples than in seawater, meaning that simple two-end-member mixing cannot adequately explain the groundwater salinization. Stable isotopes of high-TDS saline and brine samples fall between the seawater and CWS mixing lines, suggesting a potential for three-end-member mixing processes (Douglas et al., 2000). Therefore, we considered SW01 (seawater) and P18 (mostly saline, but with relatively depleted stable isotopes) as two saline end-members. P16-100, which was most likely recharged during the last deglaciation, was chosen to represent fresh end-members that could have been impacted by overlying seawater or CSW during the Holocene transgression. In Fig. 8, an inferred salinization zone is established that includes almost all saline and brine groundwater samples, demonstrating the salinization processes in which fresh groundwater mixed with either seawater, CSW, or a mixture of both.

The fresh and brackish groundwater samples, on the other hand, have low Cl⁻ concentrations and less ¹⁸O, so they deviate from the assumed salinization zone and approach the river samples in Fig. 8, implying a river water–groundwater mixing trend. LH02 (with low δ¹⁸O) and SH02 (with high δ¹⁸O) were selected to represent river water end-members for different river recharge of groundwater in the study area, while G09-15 (saline but with river-like stable isotope levels) was considered as a groundwater end-member. It is presumed that a freshening zone could form between the two river water-groundwater mixing lines, indicating the occurrence of freshening processes, which would be in agreement with the continental runoff dilution discussed in Sect. 5.2.

Figure 8 Plot showing the relationship between Cl⁻ and δ¹⁸O in different water samples as a means to probe various mixing processes in the Luanhe River delta. The symbols are the same as used in Fig. 6. The green area is the assumed freshening zone, and the purple area is the assumed salinization zone.

6 Interpretation of paleo-environmental development

Based on an analysis of a range of evidence related to Quaternary geographic evolution, it is possible to understand the changes in hydrogeological conditions in the past (van Engelen et al., 2018). Other authors (Wang et al., 1981; Peng et al., 1981; Li et al., 1982, 1983; Xu et al., 2018) found that the Pleistocene transgression events – related to marine isotope stage (MIS) 3 and 5 – reached the study area once, which would have resulted in groundwater salinization. Since the last deglaciation (about 15 ka BP), the paleo-coastline has dropped to a depth of approximately 100 m below the present sea level along the shelf edge (Li et al., 2014). Stronger river downcutting and flushing in the study region would have promoted a large fresh recharge of groundwater. For example, P16-100 (freshwater) was sampled from a relatively deep position (100 m below the surface) and has an estimated groundwater age of between 15 959 to 17 490 a BP, which likely provides evidence that the salinization groundwater related to MIS 5 and/or 3 marine transgression could have been flushed out until the latest Late Pleistocene. Accordingly, we believe that the observed saline groundwater in the Luanhe River delta is probably related to the subsequent Holocene marine transgression. This research illuminates the evolutionary pattern of saline groundwater, as shown in Table 1 and Fig. 9, where three phases are synthesized and reconstructed.

Table 1 Saline groundwater evolution processes in the study area.

The global sea level was affected by deglaciation of the ice sheet (Fairbanks, 1989), causing the sea level to rise rapidly during the deglaciation period (15.4–7 ka BP) (Li et al., 2014). The Holocene transgression stage, which occurred between 9 and 7 ka BP, resulted in the study area being inundated by seawater (Xu et al., 2015; Xue, 2009, 2014) (Fig. 9a). On the one hand, there would have been a tendency for the denser seawater to infiltrate through the aeration zone (Santucci et al., 2016); on the other hand, the sea-level rise would have caused the seawater–freshwater interface to move landward (Ferguson and Gleeson, 2012), both of which contributed to paleo-seawater intrusion. G08-40 contains 27.173 g L⁻¹ of TDS, which is similar to that of SW01. Simultaneously, the residence time (9810–6884 a BP) indicates the existence of paleo-seawater trapped by low-permeability aquitard sediments, which may be another critical salinity source for neighboring aquifers in the coastal zone (Post and Kooi, 2003; Lee et al., 2016).

Figure 9 Diagram of paleo-environmental development during the Holocene and evolutionary pattern of saline groundwater.

The presence of paleo-seawater intrusion during the Quaternary has been recorded in other coastal regions worldwide (Groen et al., 2000; Bouchaou et al., 2009; Wang and Jiao, 2012; Delsman et al., 2014; Tran et al., 2020; Han et al., 2020). According to the works described above, the salinity of the groundwater after salinization could not exceed that of seawater due to paleo-seawater intrusion.

Other salinization processes that occurred during paleo-environmental growth are likely to be correlated with this brine groundwater.

6.1 Phase 2: highstand system tract–old Luanhe River delta development (7–3.5 ka BP)

The good fit between the measured hydrochemistry and simulated evaporation lines (Figs. 6 and 7) is an indicator that the brine samples were associated with the seawater that was exposed to evaporation during geological history. Previous research has revealed that the lagoon was active during the progradation of the old Luanhe River delta between 7 and 3.5 ka BP (He et al., 2020; Xu et al., 2020). Meanwhile, the region has had a relatively arid climate since 5500 a BP, which may have led to increased evaporation (Jin, 1984). The ancient lagoon would have been an ideal location for evaporating seawater that had been trapped due to storms or tides (Fig. 9b). As a result, concentrated saline water (CSW) with a salinity greater than that of seawater would have been created, and the CSW would have undergone two processes: (1) it would have infiltrated and descended to the lower part of the aquifer due to its higher density, before combining with the salinized groundwater from phase 1, resulting in a three-end-member mixing scenario in the relationship diagram (Fig. 8); (2) after reaching saturation during the later stages of evaporation, precipitation of minerals such as gypsum, calcite, and halite would have occurred, which would have been redissolved by meteoric waters or seawater, resulting in high-salinity water which would then have been subjected to the above process. The Br $/$ Cl ratios in certain fresh or brine groundwater samples deviate from the evaporation line (Fig. 7b), which may be related to halite precipitation and redissolution. These two processes caused the groundwater salinity to rise even further, resulting in the formation of brine groundwater with 3 times the TDS of seawater, such as G03-20, with a residence time range of 4323–5590 a BP.

6.2 Phase 3: new Luanhe River delta development (3.5 ka BP to present)

Since about 3500 a BP, a nearly 90^∘ diversion of the Luanhe River channel in the study area has resulted in new delta development (Wang et al., 2007; Xue, 2016). There are some signs of a lagoon environment in the new Luanhe River delta (Cheng et al., 2020), and, as previously discussed, the brine groundwater sample G10-30 can be attributed to evaporation in a lagoon setting (Fig. 9c). However, some factors are likely to limit CSW formation in the study area: (1) the relatively low evaporation capacity due to the semi-humid climate since about 2.5 ka BP (Jin, 1984), (2) diluvial deposits or artificial reclamation would have filled the coastal lowlands such as lagoons, and (3) offshore levees prevent the seawater from flooding inland during storms or tides. Unlike the old Luanhe River delta, these factors may also explain why the current Luanhe River delta does not have high-TDS brine groundwater.

In addition, the brackish and low-TDS saline groundwater with a relatively modern age (e.g., G09-15) and the river-like stable isotope levels (Figs. 4 and 8) are compelling evidence that freshening processes have occurred in the delta plain. With the semi-humid paleoclimate, some abandoned channels have developed into small rivers such as the Suhe River and Shahe River following the diversion of the ancient Luanhe River (Gao, 1981). Firstly, the lateral recharge from the surface stream plays a role in washing out the salty groundwater. Secondly, due to the unsuitability of saline groundwater throughout human history, river irrigation has been commonly used for agricultural activities in the study region, which has freshened the upper saline aquifer (Fig. 9c). Some groundwater samples found above the seawater mixing line in the Ca-Cl and SO₄-Cl relationship diagrams (Fig. 6c, d) may be related to mineral dissolution during river water or irrigation recharge. However, saline groundwater can be washed out over time in coastal zones with low-permeability marine layers and a low hydraulic gradient (van Engelen et al., 2019; Han et al., 2020).

In summary, the evolution of saline groundwater in the study area results from paleo-environment development such as sea-level changes, paleogeography, and paleoclimate, and is significantly affected by human activities. The coastal brine groundwater found on the Bohai Sea coast, such as Bohai Bay (Li et al., 2017) and Laizhou Bay (Han et al., 2014), is a special product of geological evolution. The change in sea level during the Late Pleistocene would have favored marine intrusion and a similar sedimentary environment on the Bohai coast, allowing this study to infer the following conditions for brine formation in this region: (1) stable evaporative environments (e.g., a lagoon), (2) suitable climatic conditions (e.g., aridity), (3) entry of seawater into evaporative environments (e.g., via storms or tides), and (4) a long time period for salinity accumulation.

7 Conclusions

In this study, we used a range of isotopic–geochemical methods to analyze the groundwater recharge and salinity sources in the Luanhe River delta. The isotopic results (²H, ¹⁸O, ¹⁴C) show that deep confined groundwater was recharged during the Late Pleistocene cold period, shallow saline and brine groundwater was recharged during the warm Holocene period, and shallow brackish and fresh groundwater was mainly recharged by surface water. The hydrogeochemical modeling (PHREEQC) results show that seawater or evaporated seawater is the primary source of the salt in salinized groundwater. Variation among multiple water samples in the ¹⁸O–Cl relationship further indicates multiple-end-member mixing, which is useful for assessing the salinization and/or freshening processes in the aquifers. Our study shows that multiple water types are particularly associated with complex geographic evolution in coastal areas. Rising sea levels caused lowland coastal areas to be inundated by seawater, inducing paleo-seawater intrusion. Coastal deltas developed after a significant drop in sea levels. The concentration of saline water in the lagoon environment at the front of the delta continuously added salinity to the groundwater. Thus, brine groundwater was formed under the effects of evaporation, mixing, and dissolution. In contrast, the lateral recharge by surface water and irrigation return would have caused the slow washing out of salinized groundwater in the delta plain.

Given that most coastal zones around the world have experienced transgression/regression events in the Quaternary Period, this work's findings will promote a better understanding of the origin of salinization in coastal aquifers. In addition, it is important to recognize the potential for connate saline groundwater previously preserved in adjunct aquifers to leak due to the overextraction of deep groundwater. To effectively prevent pollution from saline groundwater movement, this study recommends extensive characterization of the groundwater interface dynamics, such as the fresh/saline, fresh/brine, and brine/seawater interfaces, and that continuous monitoring of water quality and levels across the aquifers should be maintained.