The muddy coastal zone in northern Shandong is a typical chalky-gravelly coastal plain with a unique geographic location and geologic-historical background (Chen et al., 2019; Sun et al., 2023a), and its multilayered groundwater aquifer structure and endowed with a variety of hydrochemical types indicate an ideal location for studying the groundwater environment of the coastal zone. The study area, located in the coastal plain of northern Shandong, China, bordering the Bohai Sea, represents a typical coastal region. Coastal environments in the coastal aquifer of Israel and Spain (Sivan et al., 2005; Martínez-Pérez et al., 2022) use similar groundwater sampling methods (Sheng et al., 2023). The study area encompasses the alluvial plain of the Huanghe River Delta and the coastal alluvial seafloor in Laizhou Bay, characterized by a semi-humid monsoon climate with an annual average temperature of 12.1°C, 599 mm of precipitation, and 1,400 mm of evaporation. Major rivers, including the Weihe, Yuhe, Bailang, Mihe, Xiaoqing, and Huanghe Rivers, traverse the study area, providing freshwater and coarse-grained sediment to the coastal groundwater (Sun et al., 2006).

The Xiaoqing River delineates the hydrogeological unit, separating the region into the southern coast of Laizhou Bay and the Yellow River Delta. The coastal zones on either side of the Xiaoqing River show different topography, geomorphology, depositional environments and erosion levels due to riverine and oceanic influences (Chen et al., 2002; Ren et al., 2023). The southern coast of Laizhou Bay, bordered by the Bohai Sea and the Shandong Peninsula, features a multi-source, short-river deltaic depositional system with Quaternary subsurface brines predominantly concentrated in the Mahe-Weichang River basin, exhibiting a distinctive ‘dumbbell–shaped’ distribution (Gao et al., 2015). The Huanghe River, originating from the Tibetan Plateau, has formed a fan-shaped delta with an area of about 6,000 km² (Zheng et al., 2005), with sediment accumulation in the estuary creating favorable conditions for groundwater storage (Xue, 1993).

The study area’s stratigraphy consists of four layers: Early Pleistocene (burial depth: 309~172 m), Middle Pleistocene (172~45 m), Pleistocene (45~20 m), and Holocene (<20 m) (Xin and He, 1991; Guo et al., 2016) ( Table 1 ). Since the Late Pleistocene, three sea-invasion-recession events have occurred, forming three marine stratigraphic layers: Middle Holocene Huanghua (Qh², 7–2.5 ka B.P.), Late Pleistocene Xianxian (Qp₃ ³, 40–21 ka B.P.), and early Late Pleistocene Cangzhou (Qp31, 110–70 ka B.P.) (Gao et al., 2015; Gao et al., 2016). The Quaternary stratigraphy, influenced by sea-land interactions, comprises numerous weakly permeable layers and aquifers with varying hydraulic properties (Xue et al., 2000). The brine aquifer, with a multi-layered structure ( Figure 1b ), is primarily developed in three marine strata within Quaternary loose sediments (Zheng et al., 2014; Liu et al., 2017).

This study, supported by the 2023 Joint Fund ‘Mechanism of Underground Brine Resources Formation and Evolution in Shandong’s Coastal’, collected water samples from various water bodies in the study area. Our collection totaled 152 samples, compiled from 138 groundwater samples, 12 surface river samples, a single Bohai Sea sample, and a local rainfall sample, with the sea encroachment line encompassing most sample locations ( Figure 1a ). Primarily, the groundwater sampling was accomplished by utilizing local civil wells, machine wells, and monitoring wells to gather shallow (from SHA) and deep groundwater samples (from DPA). For detailed analysis of brine variances, deep brines in Laizhou Bay were categorized into Deep-upper brine (Late Pleistocene) and Deep-lower brine (Early Pleistocene). All samples were immediately filtered on-site using a 0.22 µm membrane. Analytic cation samples were acidified to a pH below 2 with 6 mmol/L nitric acid and stored in 500 mL plastic sampling bottles for water chemistry analysis. Additionally, hydrogen and oxygen isotope testing was conducted on 141 samples using 10 mL brown glass sampling bottles.

In situ measurements of groundwater’s redox potential, pH, conductivity, and temperature were performed using a Manta+ water quality multi–parameter analyzer produced by Eureka, U.S.A. The concentrations of various ions namely K⁺, Na⁺, Ca²⁺, Mg²⁺, and Br^– were tested using an Inductively Coupled Plasma analyzer (ICAP–7400). The evaluation of SO₄ ^2– and Cl^– ions was done using Chromatographic Analysis (ICS–600), while the HCO₃ ^– ion concentration were established via titration. The Marine Geological Testing Centre of China’s Ministry of Natural Resources employed High Temperature Pyrolysis–Isotope Ratio Mass Spectrometry (HTP–IRMS) to measure the water samples’ stable isotope contents (δ²H, δ¹⁸O), using the Vienna Standard Marine Water (V–SMOW) as the calculating index. The isotopic uncertainties of δ²H and δ¹⁸O were ±1.0‰ and ±0.2‰, respectively. Additionally, Beta Laboratory performed tests and analysis on the stable isotopes (δ²H, δ¹⁸O) of the groundwater samples.

In the study region, TDS contents in SHA and DPA groundwater ranged from 0.38–34.20 g/L and 0.46–156.96 g/L, respectively, with average values of 8.06 g/L and 72.07 g/L. The dominant hydrochemical compositions were Cl–Na and Cl–Na·Mg ( Figure 2 ). Major ions indicated a transition from freshwater to brine, with cations evolving from Ca²⁺ and Mg²⁺ to Na⁺, and anions from Bicarbonate and Sulfate to Chloride. The hydrochemistry of freshwater in SHA is mainly HCO₃–Ca and Cl·SO₄–Ca·Mg, while in DPA, it’s predominantly HCO₃–Ca, HCO₃–Ca·Na. The brackish water mainly comprises of Cl·SO₄–Na and Cl·SO₄–Ca·Mg, while Cl–Na characterizes both saline and brine waters. DPA groundwater exhibited a larger range of TDS contents and more complex hydrochemical types. Most brines were similar to local seawater, suggesting a seawater origin.

From a regional perspective, Cl–Na type groundwater in the offshore DPA of the Huanghe River Delta is more homogeneous, likely due to the extensive marine layer. In contrast, groundwater along the southern coast of Laizhou Bay shows more diverse hydrochemical types, indicating potential mixing with other aquifers.

TDS contents of groundwater samples generally decreased with increasing distance from the coastline ( Figure 3 ). Within the coastal (<20 km) area, DPA groundwater had higher TDS contents, often exceeding local seawater (24.1 g/L). Beyond 20 km, SHA groundwater TDS significantly lowered to freshwater levels. Near the sea, SHA groundwater was mainly saline, with some brine levels within 10 km. In the 20~30 km range, only DPA groundwater had relatively higher TDS, similar to seawater, suggesting past seawater intrusion influence. The fitted curve showed DPA groundwater had higher TDS than SHA groundwater, with greater change near the coastline. SHA groundwater near the coast had TDS equivalent to local seawater, indicating significant seawater mixing and intrusion influence in near-coastal Holocene strata.

The relationship between ions and Cl^- can indicate the source of groundwater salinity (e.g. Edmunds et al., 2006). In ion vs. Cl^- plots ( Figure 4 ), most water samples, particularly saline and brine waters, align closely with the seawater ratio line, as seen in Na⁺ vs. Cl^- ( Figure 4a ), Mg²⁺ vs. Cl^- ( Figure 4b ), SO₄²^- vs. Cl^- ( Figure 4c ), and Br^- vs. Cl^- relationships ( Figure 5a ). Fresh and brackish water samples, however, are scattered and shifted upward. Lower salinity groundwater is likely influenced by mineral weathering, which dissolves certain ions and increases total solutes. In high-salinity conditions, re-dissolving mineral ions is difficult. Using the mass balance model (Equation 2), the mixing line between end-element samples was calculated and plotted across water chemistry indicators (Figure 7 ). Most ions align with seawater and brine end-members, but Ca²⁺ vs. Cl^- shows a scattered distribution for groundwater samples ( Figure 4d ), with significant blending of SHA and DPA freshwater end-members.

Saline water with TDS contents exceeding seawater and brine is closely distributed near the seawater evaporation line, indicating groundwater salinity may be linked to seawater or concentrated saline water. Na⁺, Mg²⁺, and SO ₄ ²^- contents in Northern Shandong’s coastal groundwater generally correlate positively with Cl^- concentration. Cl^-, Na⁺, and Mg²⁺ are key seawater components. Groundwater salts could originate from seawater intrusion, marine aerosols, or dissolution of salt rocks by precipitation (see Figures 4a–c ). Most fresh, brackish, and saline water lies below the mixing line ( Figure 4a ), showing Na⁺ depletion and possible cation exchange. Brackish and saline waters deviating from the mixing line may exhibit SO₄²^- enrichment ( Figure 4c ), suggesting sulfate mineral dissolution. Further hydrogeochemical analyses are needed to confirm these observations and the groundwater formation evolution due to aquifer sediment contact.

Located in a temperate semi-humid monsoon climate zone, the study area is a transitional region between land and sea with high evaporation rates. Using the global atmospheric precipitation line equation (δ²H = 8δ¹⁸O + 10, Craig, 1961), the local meteorological water line (LMWL) was determined as δ²H = 7.8δ¹⁸O + 6.3, based on monthly average rainfall data in Yantai from 1986–1990. The δ¹⁸O and δ²H values for average monthly rainfall were obtained from the IAEA Global Network of Isotopes in Precipitation (GNIP) (IAEA/WMO, 2006). The regression equation for water samples had a slope of 5.8 (R² = 0.91) ( Figure 7 ), lower than global and local atmospheric precipitation lines. The Bohai Coastal Plain, affected by past climatic conditions and sea intrusion-recession events, shows distinct shifts in δ¹⁸O and δ²H values in groundwater (–78.07 ‰ to –19.17 ‰ for δ²H, –11.26 ‰ to –2.09 ‰ for δ¹⁸O). Most groundwater samples are isotopically below global and local atmospheric precipitation lines, especially those with δ¹⁸O > –8 ‰. Isotopes in saline and brine water are enriched, while fresh and brackish water samples align with river water.

Most groundwater samples experience evapotranspiration, resulting in smaller regression equation slopes and intercepts, as depicted by the green dashed line in Figure 8 . This is evident in saline and brine waters with significantly depleted δ¹⁸O and δ²H values. Fresh, brackish, and partially saline waters are near the atmospheric precipitation line, indicating minimal evapotranspiration impact. While saline/brine water from seawater evaporation typically has δ¹⁸O and δ²H values above standard seawater (δ²H = 0 ‰, δ¹⁸O = 0 ‰) (Horita, 2005), study area groundwater has lower values, suggesting alternative recharge sources. Subsurface brine samples, particularly SHA brine, are depleted in δ¹⁸O and δ²H, indicating freshwater mixing with terrestrial isotopes during brine formation. SHA groundwater samples (HD38 and HD40) near the Xiaoqing River’s offshore zone have high δ¹⁸O and δ²H values, suggesting post-evaporation seawater origin within SHA, similar to standard seawater. Climate factors cause isotopic depletion in most DPA brines. Notably, δ¹⁸O and δ²H in the Huanghe River Delta DPA brine sample (HD46, ~260 m) closely match the Huanghe River water sample, indicating a strong link between river water and Early Pleistocene strata brine moisture sources.

As previously discussed, groundwater in the study area results from the mixing of seawater, freshwater, and brine. A three-end-member mixing model (SHA freshwater (LB01) and DPA freshwater (LB32) as freshwater end–members, seawater (Sea) as the seawater end–member, and DPA brine (LB44) as the brine end–member) using Cl^– concentration and δ¹⁸O values was applied to trace mixing trajectories and quantify variations in water types ( Figure 7 ). The model shows a correlation between Cl^– concentration and δ¹⁸O values, with δ¹⁸O increasing as Cl^– concentration rises, indicating the model’s validity. Most samples are above the DPA fresh-brine mixing line, with freshwater samples between the two mixing lines, having δ¹⁸O values from –9.5 to –8.2 ‰. The mass balance model (Equation 2) can estimate the contributions of freshwater, seawater, and brine to the mixed samples.

River water from the Bailang River (BLH) and Weihe River (WH) has high δ¹⁸O values due to upstream reservoir mixing and evaporation. SHA groundwater, recharged by these rivers, shows enriched δ¹⁸O values near the river basins. Post-sea recession interaction between rivers and groundwater is a key process for flushing saline water and reducing groundwater salinity (Kwong and Jiao, 2016). Atmospheric precipitation also contributes to high δ¹⁸O in SHA water. Most Huanghe River Delta groundwater is in the mixing zone. SHA brines (HD38, HD40) near the Xiaoqing River have high Cl^– –δ¹⁸O values, suggesting a blend of evaporated seawater or modern precipitation. This differs from the Laizhou Bay south coast groundwater origin, hypothesizing brine from Holocene seawater via strong evaporation. DPA brine (HD46) at 260m depth has δ¹⁸O values similar to the Huanghe River, indicating Early Pleistocene brine recharge by river water. Other DPA brines at 100m depth mix with desalinated groundwater and seawater, clustering within the mixing line.

In Section 4.2, we discussed the water–rock interactions of groundwater. To further understand these interactions, we conducted inversion simulations to identify hydrochemical reactions and changes in ion concentrations.

We selected three paths along the groundwater flow direction to simulate the south bank of Laizhou Bay and the Huanghe River Delta. The paths in the Laizhou Bay area are Path I (LB30 → LB35), Path II (LB35 → LB73), and Path III (LB73 → LB46). In the Huanghe River Delta, we considered Path IV (HD23 → HD59) and Path V (HD59 → HD62) ( Table 2 ). Given the selected pathways and the hydrogeological background of the study area, we finally selected dolomite, calcite, gypsum, rock salt, quartz, sodium feldspar, potassium feldspar, and cation exchange as the ‘probable mineral phases’ in the Huanghe River delta (Tian et al., 2021), with the additional presence of hematite and pyrite in the Laizhou Bay area (Wang et al., 2023; Li et al., 2021) ( Table 3 ). Dolomite was set to only precipitate, while gypsum and calcite were set to only dissolve. Simulations accounted for cation exchange between Ca²⁺, Mg²⁺, and Na⁺ under seawater intrusion or aquifer renewal conditions. The outcomes were used to analyze CO₂ dissolution and release, identifying inorganic carbon sources and sinks. The analysis focused on seven elements: Na⁺, K⁺, Ca²⁺, Mg²⁺, Cl^–, HCO₃ ^–, and SO₄ ^2–.

The cation exchange involving Na⁺ replacing Ca²⁺ or Mg²⁺, and Ca²⁺ replacing Mg²⁺ in all pathways along the south coast of Laizhou Bay aligns with the findings in Section 4.2.2, suggesting active seawater intrusion. In Path I, the transition from freshwater to saltwater involves the dissolution of gypsum, rock salt, and potassium salt, leading to increased ion content, dolomite precipitation, and CO₂ release. Despite Ca²⁺ consumption by dolomite precipitation, Ca²⁺ concentration in groundwater slightly rises due to gypsum dissolution and Na⁺ exchange. Path II (from saltwater to brine) and Path III (from brine to brine) exhibit similar patterns, with calcite, gypsum, rock salt, and potassium salt dissolution and CO₂ absorption. High Na⁺ concentration in coastal groundwater triggers Na⁺ replacement reactions with Ca²⁺ and Mg²⁺ in the aquifer. Extensive gypsum dissolution increases calcium ions, enhancing calcite solubility and leading to calcite dissolution.

The Huanghe River Delta includes Path IV (from brackish water to brine) and Path V (from high concentration (TDS content: 121 g/L) brine to low concentration (TDS content: 92 g/L) brine). Path IV involves gypsum dissolution and potassium salt precipitation. Path V involves the dissolution of albite and the precipitation of potassium feldspar, with CO₂ absorption. Seawater intrusion at the HD62 sample site leads to albite dissolution and potassium feldspar precipitation, indicating significant interaction between seawater and subsurface brine. Both paths show equal cation exchange between Ca²⁺ and Mg²⁺, due to the ionic radius and adsorption energy of Ca²⁺ in the coastal aquifer, facilitating exchange with clay minerals (Capuano and Jones, 2020). In the saline environment, Na⁺ and Mg²⁺ concentrations are elevated, with Na⁺ being predominant (Qin et al., 2023).

Significant fluctuations in the Br^– vs. Cl^– ratio in freshwater and brackish water with lower Cl^– concentration are noted, influenced by atmospheric precipitation or anthropogenic factors ( Figure 5b ). Offshore groundwater samples are close to the seawater ratio line, indicating predominant seawater mixing as the source of salinity. The Br^–/Cl^– ratio is unreliable for identifying salinity sources when TDS is less than 2 g/L (Han et al., 2011). SHA brine and certain brackish waters show lower Br^–/Cl^– values compared to seawater, enclosed within the black dashed line of the study region, which is due to the significant increase in Cl^– concentration from the dissolution of rock salt. Highly mineralized SHA brine and Deep–upper brine (LZ) show a linear decrease in Br^–/Cl^– with increasing Cl^–.

Typically, the Br^–/Cl^– ratio stays consistent during seawater evaporation, that is, as groundwater salinity intensifies (Bottomley et al., 1994; Vengosh and Hendry, 2001), until it reaches rock salt saturation. Underground brine samples from the Huanghe River Delta, particularly the DPA brine, exhibit larger Br^–/Cl^– ratios, exceeding those of the seawater ratio line. This suggests that they may have passed the rock salt saturation stage, advancing towards the lagoonal salt saturation phase, which could have prompted Cl^– precipitation in water, followed by mixing with freshwater, subsequently reducing salinity. Deep–lower brine (LZ) samples within the red dotted line, located along the river–adjacent coastal, with a burial depth of 80~200m, possess higher Br^–/Cl^–. This could be attributed to the decomposition of organic matter within the marine layer of the DPA groundwater.

The HFE–Diagram analysis confirms that groundwater salinity is high (within a 20 km radius from the coastline), with significant seawater intrusion. The mineralization of SHA and DPA groundwater decreases with desalination and increases with seawater intrusion in Figure 12 . Brine samples cluster near the end stages of intrusion. The desalination endpoints in SHA align with LB01 ( Figure 12a ), and in DPA with LB32 ( Figure 12b ), consistent with previous mixing process findings in Figure 6 . In SHA, 63.77% of samples undergo desalination and 36.23% intrusion, with Na–Cl as the predominant phase (55.07%). This reflects the influence of shallow water controlled by hydrologic cycle processes and surface water recharge. In DPA, 79.31% of samples experience intrusion and 20.69% desalination, with Na–Cl being most abundant (89.66%). The DPA groundwater, on the other hand, is predominantly saline, which is controlled by tectonic conditions of confinement and rock salt dissolution. Extensive groundwater extraction has impacted both quantity and quality (Currell et al., 2012), but most samples show a trend towards seawater intrusion, reflecting global coastal aquifer trends (Ri et al., 2022; Azari and Tabari, 2024).