Delineating multiple salinization processes in a coastal plain aquifer, northern China: hydrochemical and isotopic evidence

Groundwater is an important water resource for agricultural irrigation and urban and industrial utilization in the coastal regions of northern China. In the past 5 decades, coastal groundwater salinization in the Yang–Dai river plain has become increasingly serious under the influence of anthropogenic activities and climatic change. It is pivotal for the scientific management of coastal water resources to accurately understand groundwater salinization processes and their causative factors. Hydrochemical (major ion and trace element) and stable isotopic (δ¹⁸O and δ²H) analysis of different water bodies (surface water, groundwater, geothermal water and seawater) were conducted to improve understanding of groundwater salinization processes in the plain's Quaternary aquifer. Saltwater intrusion due to intensive groundwater pumping is a major process, either by vertical infiltration along riverbeds which convey saline surface water inland, and/or direct subsurface lateral inflow. Trends in salinity with depth indicate that the former may be more important than previously assumed. The proportion of seawater in groundwater is estimated to have reached up to 13 % in shallow groundwater of a local well field. End-member mixing calculations also indicate that the geothermal water with high total dissolved solids (up to 10.6 g L⁻¹) with depleted stable isotope compositions and elevated strontium concentrations (> 10 mg L⁻¹) also mixes locally with water in the overlying Quaternary aquifers. This is particularly evident in samples with elevated Sr ∕ Cl ratios (> 0.005 mass ratio). Deterioration of groundwater quality by salinization is also clearly exacerbated by anthropogenic pollution. Nitrate contamination via intrusion of heavily polluted marine water is evident locally (e.g., in the Zaoyuan well field); however, more widespread nitrate contamination due to other local sources such as fertilizers and/or domestic wastewater is evident on the basis of NO₃ ∕ Cl ratios. This study provides an example of how multiple geochemical indicators can delineate different salinization processes and guide future water management practices in a densely populated water-stressed coastal region.

1 Introduction

Coastal regions are key areas for the world's social and economic development. Approximately 40 % of the world's population lives within 100 km of the coast (UN Atlas, 2010). Worldwide, these areas have become increasingly urbanized, with 14 of the world's 17 largest cities located along coasts (Creel, 2003). China has 18 000 km of continental coastline, and around 164 million people (approximately 12 % of the total population) living in 14 coastal provinces; nearly 80 % of these people inhabit three coastal “economic zones”, namely Beijing–Tianjin–Hebei, the Yangtze River delta and the Pearl River delta (Shi, 2012). The rapid economic development and growing population in these regions have greatly increased demand for fresh water. Meanwhile, they are also confronted with increased sewage and other wastewater discharge into coastal environments.

Groundwater resources play crucial roles in the social, economic and ecologic function of global coastal systems (IPCC, 2007). Coastal aquifers connect with the ocean and with continental hydroecological systems (Moore, 1996; Ferguson and Gleeson, 2012). As an important freshwater resource, groundwater may be overextracted during periods of high demand, which are often periods of low recharge and/or surface water availability (Post, 2005). Overexploitation of groundwater can therefore result in seawater intrusion, as well as related environmental issues such as land subsidence. Seawater intrusion has become a global issue and related studies can be found from coastal aquifers around the world, including Israel (Sivan et al., 2005; Yechieli et al., 2009; Mazi et al., 2014), Spain (Price and Herman, 1991; Pulido-Leboeuf, 2004; Garing et al., 2013), France (Barbecot et al., 2000; de Montety et al., 2008), Italy (Giambastiani et al., 2007; Ghiglieri et al., 2012), Morocco (Bouchaou et al., 2008; El Yaouti et al., 2009), the USA (Gingerich and Voss, 2002; Masterson, 2004; Langevin et al., 2010), Australia (Zhang et al., 2004; Narayan et al., 2007; Werner, 2010), China (Xue et al., 2000; Han et al., 2011, 2015), Vietnam (An et al., 2014), Indonesia (Rahmawati et al., 2013), India (Radhakrishna, 2001; Bobba, 2002) and Brazil (Montenegro et al., 2006; Cary et al., 2015). Werner et al. (2013) provides a comprehensive review of seawater intrusion processes, investigation and management.

Seawater, or saltwater, intrusion is a complicated hydrogeological process, due to the impact of aquifer properties, anthropogenic activities (e.g., intensive groundwater pumping, irrigation practices), recharge rates, variable density flow, tidal activity and effects relating to global climate change, such as sea level rise (Ghassemi et al., 1993; Robinson et al., 1998; Smith and Turner, 2001; Simpson and Clement, 2004; Narayan et al., 2007; Werner and Simmons, 2009; Wang et al., 2015). Understanding the complex interactions between groundwater, surface water and seawater is thus essential for effective management of coastal water resources (Mondal et al., 2010). Very different salinization patterns may arise as a result of diverse interactions in coastal settings (Sherif and Singh, 1999; Bobba, 2002; Westbrook et al., 2005). Modeling has shown that generally, seawater intrusion is more sensitive to groundwater pumping and recharge rates in comparison to tidal fluctuation and sea level rise (Narayan et al., 2007; Ferguson and Gleeson, 2012). However, most models of seawater intrusion require simplification of the coastal interface zone. Relatively few studies have focused on delineating complex interactions among the surface–ground–seawater continuum in estuarine environments, including the effects of vertical infiltration of seawater into aquifers through river channels, compared to subsurface lateral landward migration of the freshwater–saltwater interface. Recent data indicate that such processes may be more important in causing historical salinization of coastal groundwater than previously appreciated (e.g., Cary et al., 2015; Lee et al., 2016; Larsen et al., 2017).

Additionally, groundwater in coastal aquifers may be affected by other salinization processes, such as input of anthropogenic contaminants or induced mixing with saline water from deeper or adjacent formations, which may include the high total dissolved solids (TDSs) in geothermal water or brines emplaced in the coastal zone over geologic history. The data from China's marine environment bulletin released on March 2015 by the State Oceanic Administration showed that the major bays, including Bohai Bay, Liaodong Bay and Hangzhou Bay, are seriously polluted, with inorganic nitrogen and active phosphate being the major pollutants (SOA, 2015). Seawater intrusion in China is most serious in the circum-Bohai Sea region (Han et al., 2011, 2016a); and due to the heavy marine pollution, the impacts of anthropogenic activities on groundwater quality in future may not simply be a case of simple saltwater intrusion. This region is also characterized by deep brines and geothermal waters (e.g., Han et al., 2014), which may migrate and mix with fresher groundwater due to intensive water extraction. Depending on the specific processes involved, additional contaminants may mix with fresh groundwater resources in parallel with seawater intrusion, and it is thus likely to be more difficult to mitigate and remediate groundwater pollution.

A variety of approaches can be used to investigate and differentiate seawater intrusion and other salinization processes, including time-series water level and salinity measurements, geophysical methods, conceptual and mathematical modeling as well as geochemical methods (see reviews by Jones et al., 1999; Werner et al., 2013). Geochemical techniques are particularly valuable in areas where the dynamics of saline intrusion are complicated and may involve long-term processes predating accurate water level records, or where multiple salinization processes may be occurring simultaneously. These techniques typically employ the use of major ion ratios such as Cl ∕ Br and Cl ∕ Na, which are indicative of solute origins (Edmunds, 1996; Jones et al., 1999). Other ionic ratios, involving Mg, Ca, Na, HCO₃ and SO₄, and characterization of water “types” can also be useful in determining the geochemical evolution of coastal groundwater, for example, indicating freshening or salinization, due to commonly associated ion exchange and redox reactions (Anderson et al., 2005; Walraevens, 2007). Trace elements such as strontium, lithium and boron can provide additional valuable information about sources of salinity and mixing between various end-members, as particular waters can have distinctive concentrations (and/or isotopic compositions) of these elements (e.g., Vengosh et al., 1999). Stable isotopes of water (δ¹⁸O and δ²H) are also commonly used in such studies, as they are sensitive indicators of water and salinity sources, allowing seawater to be distinguished from other salt sources (e.g., Currell et al., 2015).

This study examines the Yang–Dai coastal river plain in Qinhuangdao City, Hebei province, northern China, specifically focusing on salinization of fresh groundwater caused by groundwater exploitation in the Zaoyuan well field and surrounding areas. The study investigates groundwater salinization processes and interactions among surface water, seawater and geothermal groundwater in a dynamic environment, with significant pressure on water resources. Qinhuangdao is an important port and tourist city of northern China. In the past 30 years, many studies have investigated seawater intrusion and its influencing factors in the region using hydrochemical analysis (Xu, 1986; Yang et al., 1994, 2008; Chen and Ma, 2002; Sun and Yang, 2007; Zhang, 2012) and numerical simulations (Han, 1990; Bao, 2005; Zuo, 2009). However, these studies have yet to provide clear resolution of the different mechanisms contributing to salinization, and they have typically ignored the role of anthropogenic pollution and groundwater–surface water interaction. This study is thus a continuation of previous investigations of the region, using a range of hydrochemical and stable isotopic data to delineate the major processes responsible for increasing groundwater salinity, including lateral subsurface seawater intrusion, vertical leakage of marine-influenced surface water, induced mixing of saline geothermal water and anthropogenic pollution. The goal is to obtain a more robust conceptual model of the interconnections between the various water sources under the impact of groundwater exploitation. The results provide significant new information to assist water resources management in the coastal plain of Bohai Bay, and other similar coastal areas globally.

2 Study area

The Yang–Dai coastal river plain (Fig. 1) covers approximately 200 km² of the west side of Beidaihe District of Qinhuangdao City, northeastern Hebei Province. It is surrounded by the Yanshan Mountains to the north and west, and the southern boundary of the study area is the Bohai Sea. The plain declines in topographic elevation (with an average slope of 0.008) from approximately 390 m a.s.l. (meters above sea level) in the northwest to 1–25 m in the southeast, forming a fan-shaped distribution of incised piedmont alluvial plain sediments. Zaoyuan well field, located in the southern edge of the alluvial fan, approximately 4.3 km from the Yang River estuary, was built in 1959 (Xu, 1986) as a major water supply for the region (Fig. 1).

Figure 1 Location map (a) for showing geological background and water sampling sites in the study area, and (b) hydrogeological cross section of Yang–Dai river Plain ($P-P^{\prime }$ in a) (modified from Han, 1988). The surface area covered from the > 250 mg Cl L⁻¹ contour line to the coastline refers the seawater zones.

2.1 Climate and hydrology

The study area is in a warm and semihumid monsoon climate. On the basis of a 56-year record in Qinhuangdao, the mean annual rainfall is approximately 640 mm, the average annual temperature is approximately 11 ^∘C, and mean potential evaporation is 1469 mm. A total of 75 % of the total annual rainfall falls in July–September (Zuo, 2006), during the East Asian Summer Monsoon. The average annual tide level is 0.86 m (meters above Yellow Sea base level), while the high and low tides are approximately 2.48 and −1.43 m.

The Yanghe River and Daihe River, originating from the Yanshan Mountains, are the major surface water bodies in the area, flowing southward into the Bohai Sea (Fig. 1). The Yang River is approximately 100 km long with a catchment area of 1029 km² and average annual runoff of 1.11 × 10⁸m³ a⁻¹ (Han, 1988). Dai River has a length of 35 km and catchment area of 290 km², with annual runoff of 0.27 × 10⁸m³ a⁻¹. The rivers become full during intense rain events, and revert to minimal flow during the dry season – in part this is related to impoundment of flow in upstream reservoirs.

2.2 Geological and hydrogeological setting

Groundwater in the area includes water in Quaternary porous sediment as well as fractured bedrock in the northern platform area. Fractured rock groundwater volume mainly depends on the degree of weathering and the nature and regularity of fault zones (Fig. 1). The strata outcropping in the west, north and eastern edge of the plain include Archean, Proterozoic and Jurassic metamorphic and igneous rocks, which also underlie the Quaternary sediments of the plain (from which most samples in this study were collected). The basement faults under Quaternary cover are mainly NE-trending and NW-trending (Fig. 1); these structures control the development and thickness of the overlying sediments, as well as the distribution of hot springs and geothermal anomalies. Fault zones are thought to be the main channel for transport of thermal water from deeper to shallower depths.

The Quaternary sediments are widely distributed in the area, with the thickness ranging from approximately 5–80 m (mostly 20–40 m) up to more than 100 m immediately adjacent to the coastline. The bottom of the Holocene (Q₄) unit in most areas consists of clay, making the groundwater in the coastal zone confined or semiconfined, although there are no regional, continuous aquitards between several layers of aquifer-forming sediments (Fig. 1b). The aquifer is mainly composed of medium sand, coarse sand and gravel layers with a water table depth of 1–4 m in the phreatic aquifer, and deeper semiconfined groundwater (where present and hydraulically separated from the phreatic aquifer) hosted in similar deposits with a potentiometric surface 1–5 m below topographic elevation (Zuo, 2006).

The general flow direction of groundwater is from northwest to south, according to the topography. The main sources of recharge are from infiltration of rainfall, river water and irrigation return flow, as well as lateral subsurface inflow from the piedmont area. Naturally, groundwater discharges into the rivers and the Bohai Sea. Apart from phreatic water evaporation, groundwater pumping for agricultural, industrial and domestic usage (including seasonal tourism) are currently the main pathways of groundwater discharge.

Geothermal water discharges into shallow Quaternary sediments near the fault zones, evident as geothermal anomalies (Hui, 2009). The temperature of thermal water ranges from 27 to 57 ^∘C in this low-to-medium temperature geothermal field (Zeng, 1991). Deeper thermal water is stored in the Archaeozoic granite and metamorphic rocks; major fracture zones provide pathways into the overlying Quaternary sediments (Pan, 1990; Shen et al., 1993; Yang, 2011).

2.3 Groundwater usage and seawater intrusion history

Shallow groundwater pumped from the Quaternary aquifer occupies 94 % of total groundwater exploitation, and is used for agricultural irrigation (52 % of total groundwater use), industrial (32 %) and domestic water (16 %) (Meng, 2004). Many large and medium-sized reservoirs were built in the 1960s and 1970s, meaning that the surface water was intercepted and downstream runoff dropped sharply, even causing rivers to dry up in drought years. With the intensification of human socioeconomic activities and growing urbanization, coupled with extended drought years (severe drought during 1976–1989 in north China) (Wilhite,1993; Han et al., 2015), increased groundwater exploitation to meet the ever-growing fresh water demand resulted in groundwater table declines and seawater intrusion (SWI) in the aquifers.

The pumping rate in the Zaoyuan well field gradually increased from 1.25 million m³ a⁻¹ in the early 1960s to 3.5 million m³ a⁻¹ in the late 1970s, and beyond 10 million m³ a⁻¹ in the 1980s. During 1966–1989, the planting of paddy fields became common, resulting in significant agricultural water consumption. This caused formation of a cone of depression in the Quaternary aquifer system. Groundwater pumping in this region mainly occurs in spring and early summer, with typical pumping rates of 70 000–80 000 m³ d⁻¹. Pumping from the Zaoyuan well field occurs in wells approximately 15 to 20 m deep. Groundwater tables decline sharply and reach their lowest level during May, before the summer rains begin, and recover to their yearly high in January–February (Fig. 2). In May 1986, the groundwater table in the depression center, which is located in Zaoyuan–Jiangying (Fig. S1 in the Supplement), decreased to −2 m a.s.l. and the area with groundwater tables below sea level covered 28.2 km². The local government commenced reduction in groundwater exploitation in this area after 1992, and groundwater tables began to decrease more slowly after 1995, even showing recovery in some wells. However, during an extreme drought year (1999), increased water demand resulted in renewed groundwater table declines in the region (Fig. 2). Since 2000, the groundwater tables have responded seasonally to water demand peaks and recharge (Figs. 2 and S1).

From 1990, the rapid development of township enterprises (mainly paper mills), also began to cause groundwater overexploitation in the western area of the plain. The groundwater pumping rate for paper mills reached 55 000 m³ d⁻¹ in 2002, resulting in groundwater table depressions around Liushouying and Fangezhuang (Fig. 1). The groundwater table in the western depression associated with this pumping reached −11.6 in 1991 and −17.4 m a.s.l. in 2002. After the implementation of “Transferring Qing River water to Qinhuangdao” project in 1992, the intensity of groundwater pumping generally reduced, and the groundwater table in the depression center recovered to −4.3 m a.s.l. in July 2006.

Overall, the depression area (groundwater tables below mean sea level) was recorded as 132.3 km² in May 2004 and the shape of the depression has generally been elliptical with the major axis aligned east–west. In addition to groundwater overexploitation, climate change-induced recharge reduction has also likely contributed to groundwater table declines and hence seawater intrusion (Fig. S2). The annual average rainfall declined from 639.7 between 1954 and 1979 to 594.2 mm between 1980 and 2010, a significant decrease over the last 30 years (Zhang, 2012). As indicated in Fig. S2, the severity of seawater intrusion (indicated by changes in Cl concentration, and the total area impacted by SWI, as defined by the 250 mg L⁻¹ Cl contour) correlates with periods of below-average rainfall – indicated by monthly cumulative rainfall departure (CRD, Weber and Stewart, 2004).

Groundwater quality of the area gradually became more saline from the early 1980s, with chloride concentrations increasing year by year. As early as 1979, seawater intrusion was recorded in the Zaoyuan well field. The intrusion area with groundwater chloride concentration greater than 250 mg L⁻¹ was 21.8 in 1984, 32.4 in 1991, 52.6 in 2004 and 57.3 km² in 2007 (Zuo, 2006; Zang et al., 2010). The chloride concentration of groundwater pumped from a monitored well field well (depth of 18 m, G10 in Fig. 1) changed from 90 in 1963 to 218 mg L⁻¹ in 1978, 567 in 1986, 459 in 1995 and 1367 mg L⁻¹ in 2002 (Zuo, 2006), reducing to 812 mg L⁻¹ in July 2007. The distance of estimated seawater intrusion into the inland area from the coastline had reached 6.5 in 1991 and 8.75 km in 2008 (Zang et al., 2010). In the early 1990s, 16 of 21 pumping wells in the well field were abandoned due to the saline water quality (Liang et al., 2010). Additionally, 370 of 520 pumping wells were abandoned in the wider Yang–Dai coastal river plain during 1982–1991 (Zuo, 2006).

3 Methods

In total, 80 water samples were collected from the Yang–Dai coastal river plain, including 58 groundwater samples, 19 river water samples (from 12 sites) and 3 seawater samples, during 3 sampling campaigns (June 2008, September 2009 and August 2010). Groundwater samples were pumped from 28 production wells with depths between 6 and 110 m, including 7 deep wells with depths greater than 60 m (Fig. 1). While sampling for geochemical parameters would ideally be conducted on monitoring wells, due to an absence of these, production wells were utilized. In most cases, the screened interval of these wells encompasses aquifer thicknesses of approximately 5 to 15 m above the depths indicated in Table 1.

Table 1 Physical, hydrochemical and isotopic data of the water samples collected from the Yang–Dai coastal river plain.

Figure 2 Distribution of precipitation and dynamics of the water table in the study area. Locations of the monitoring wells (W1, W2 and W3) can be seen from Fig. 1.

In this study, sampling focused predominantly on low temperature groundwater; however, geothermal water from around Danihe was also considered a potentially important ongoing source of groundwater salinity. As such, while geothermal water samples were not accessible during our sampling campaigns (as the area is now protected), data reported by Zeng (1991) were compiled and analyzed in conjunction with the sampled wells.

Measurements of physicochemical parameters (pH, temperature and electrical conductivity (EC)) were conducted in situ using a portable meter (WTW Multi 3500i). All water samples were filtered with 0.45 µm membrane filters before analysis of hydrochemical composition. Two aliquots in polyethylene 100 mL bottles at each site were collected for major cation and anion analysis, respectively. Samples for cation analysis (Na⁺, K⁺, Mg²⁺ and Ca²⁺) were treated with 6N HNO₃ to prevent precipitation. Water samples were sealed and stored at 4 ^∘C until analysis. Bicarbonate was determined by titration within 12 hr of sampling. Concentrations of cations and some trace elements (B, Sr, Li) were analyzed by inductively coupled plasma-optical emission spectrometry (ICP-OES) in the chemical laboratory of the Institute of Geographic Sciences and Natural Resources Research (IGSNRR), Chinese Academy Sciences (CAS). Only the Sr data are reported here, as the other trace elements were not relevant to the interpretations discussed (Table 1). The detection limits for analysis of Na⁺, K⁺, Mg²⁺ and Ca²⁺ are 0.03, 0.05, 0.009 and 0.02 mg L⁻¹, respectively. Concentrations of major anions (i.e., Cl⁻, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and F⁻) were analyzed using a high-performance ion chromatograph (SHIMADZU, LC-10ADvp) at the IGSNRR, CAS. The detection limits for analysis of Cl⁻, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and F⁻ are 0.007, 0.018, 0.016 and 0.006 mg L⁻¹, respectively. The testing precision of the cation and anion analysis is 0.1–5.0 %. Charge balance errors were less than 8 %. Stable isotopes (δ¹⁸O and δ²H) of water samples were measured using a Finnigan MAT 253 mass spectrometer after on-line pyrolysis with a Thermo Finnigan TC/EA in the Stable Isotope Laboratory of the IGSNRR, CAS. The results are expressed in per mill (‰) relative to international standards (V-SMOW, Vienna Standard Mean Ocean Water) and resulting δ¹⁸O and δ²H values are shown in Table 1. The analytical precision for δ²H is ±2 ‰ and for δ¹⁸O is ±0.5 ‰ . All hydrochemical, physicochemical and isotope data are reported in Table 1.

Mixing calculations were also conducted on the basis of Cl⁻ concentrations of the samples under a conservative freshwater–seawater mixing system (Fidelibus et al., 1993; Appelo, 1994). The seawater contribution for each sample is expressed as a fraction of seawater (f_sw), using the formula according to Appelo and Postma (2005):

$$\begin{array}{}\text{(1)}& f_{\mathrm{sw}}={\displaystyle \frac{{\mathrm{C}}_{\mathrm{Cl},\mathrm{sam}}-{\mathrm{C}}_{\mathrm{Cl},\mathrm{f}}}{{\mathrm{C}}_{\mathrm{Cl},\mathrm{sw}}-{\mathrm{C}}_{\mathrm{Cl},\mathrm{f}}}},\end{array}$$

where C_Cl,sam, C_Cl,f and C_Cl,sw refer to the Cl⁻ concentration in the sample, freshwater and seawater, respectively.

4 Results

4.1 Water stable isotopes (δ²H and δ¹⁸O)

The local meteoric water line (LMWL, δ²H = 6.6 δ¹⁸O + 0.3, n = 64, r² = 0.88) is based on δ²H and δ¹⁸O mean monthly rainfall values between 1985 and 2003 from Tianjin station some 120 km southwest of Qinhuangdao City (IAEA/WMO, 2006). Due to similar climate and position relative to the coast, this can be regarded as representative of the study area. Surface water samples collected from Yang River and Dai River (n = 19) have δ¹⁸O and δ²H values ranging from −10.1 to −0.6 ‰ (mean = −5.4 ‰) and −71 to −11 ‰ (mean = −43 ‰), respectively. Stable isotope compositions for surface water appear to exhibit significant seasonal variation (Fig. S3): for Yang River samples from the relatively dry season (June 2008, n = 3) had mean δ¹⁸O and δ²H values of −3.0 and −31 ‰, respectively; samples from the wet season (August 2009 and September 2010, n = 6) had mean δ¹⁸O and δ²H values of −6.6 and −48 ‰, respectively. Dai River samples showed similar results: the dry season mean δ¹⁸O and δ²H values (n = 3) were −2.6 and −32 ‰, respectively; wet season samples (n = 7) had mean δ¹⁸O and δ²H values of −6.6 and −49 ‰, respectively (Fig. 3).

The 56 groundwater samples were characterized by δ¹⁸O and δ²H values ranging from −11.0 to −4.2 ‰ (mean = −6.5 ‰) and −76 to −39 ‰ (mean = −50 ‰), respectively. Among these, shallow and deep groundwater samples showed similar mean values, although deep groundwater samples (n = 13) showed relatively narrow overall ranges (−7.8 to −5.1 ‰ and mean = −6.3 ‰ for δ¹⁸O; −58 to −43 ‰ and mean = −50 ‰ for δ²H; Fig. 3). Slight seasonal variation was evident in the groundwater isotope compositions; shallow groundwater from the dry season (n = 12) showed δ¹⁸O and δ²H values from −7.2 to −4.2 ‰ (mean = −5.7 ‰) and δ²H values from −56 to −39 ‰ (mean = −48 ‰); while during the wet season (n = 31) δ¹⁸O and δ²H values ranged from −11.0 to −5.3 ‰ (mean = −6.9 ‰) and −76 to −43 ‰ (mean = −51 ‰), respectively. Some variability was also evident in deep groundwater compositions, although only three deep samples were collected during the dry season.

Figure 3 Stable isotope compositions of different water samples collected from the study area. LMWL: local meteoric water line; GMWL: global meteoric water line (Craig, 1961); GWL: groundwater line; SWL: surface (river) water line.

From Fig. 3, it can be seen that surface water exhibits a much wider range of δ¹⁸O and δ²H values relative to groundwater, with shallow groundwater in turn more spatially variable than deep groundwater. Water samples collected in the wet season showed wider ranges of δ¹⁸O and δ²H values relative to the dry season. Most water samples of all types plot to the right of (below) the LMWL, with some surface water samples showing similar compositions to the local seawater (Fig. 3). The local seawater plots below (more negative) than typically assumed values (e.g., VSMOW = 0 ‰) for both δ²H and δ¹⁸O, and this water appears to represent an end-member involved in mixing with meteoric-derived waters in both ground and surface water (Fig. 3).

4.2 Water salinity and dissolved ions

TDS concentrations of surface water samples from Dai River range from 0.3 to 31.4 g L⁻¹, with Na⁺ and Ca²⁺ comprising 22–78 and 4–56 % of total cations and Cl⁻ comprising 36–91 % of total anions. The composition changes from Ca•Na•Mg-Cl•HCO₃ to Na-Cl water type from upstream to downstream locations along with increasing salinity; Cl⁻ concentrations vary from approximately 70 upstream to 16 700 mg L⁻¹ near the coastline, due to marine influence. Similar variation occurs along the Yang River, where samples had TDS concentrations between 0.3 and 26.1 g L⁻¹ with increasing concentrations and proportions of Cl⁻ (63.2–14 953.5 mg L⁻¹) from upstream to downstream locations. Nitrate concentrations also range from 2.8 to 65.2 mg L⁻¹ in the surface water samples, increasing downstream.

In the early 1960s, groundwater pumped from the Zaoyuan well field exhibited Ca-HCO₃ water type and chloride concentrations of 90–130 mg L⁻¹; this was followed by rapid salinization since the 1980s (see Sect. 2.3). In this study, shallow groundwater is characterized by TDS concentrations of 0.4–4.8 g L⁻¹, with Cl⁻ (34–77 %), Na⁺ (12–85 %) and Ca²⁺ (5–69 %) being the predominant major anion and cations, respectively. Groundwater hydrochemical types vary from Ca-HCO₃•Cl, Ca•Na-Cl, Na•Ca-Cl to Na-Cl (Fig. 4). Deep groundwater has TDS concentrations between 0.3 and 2.8 g L⁻¹, dominated by Ca (up to 77 % of major cations) in the upstream area and Na (up to 85 % or major cations) near the coast, with water type evolving from Ca-Cl•HCO₃ to Ca•Na-Cl and Na•Mg-Cl (Fig. 4). At present, the TDS of groundwater from the well field reaches 3.31 g L⁻¹ with Na-Cl water type (see well G15). The highest observed mixing proportions of seawater occur in shallow well G10 and deep well G2, respectively, with calculated f_sw values (according to Eq. 1) of 12.95 and 5.35 %, respectively.

Figure 4 Piper plot of different water samples. Groundwater: gw.

Hydrochemical features of thermal water from the Danihe–Luwangzhuang area (Fig. 1a) are distinct from the normal-to-low temperature groundwater. Previous work by Zeng (1991) and Hui (2009) identified geothermal water with high TDS in the fractures of deep metamorphic rock. The geothermal water was characterized by TDS values between 6.2 and 10.6 g L⁻¹ and Ca•Na-Cl water type, while Cl⁻ concentrations ranged from 5.4 to 6.5 g L⁻¹ and Sr concentrations from 6.73 to 89.8 mg L⁻¹. Some normal-to-low temperature groundwater samples collected in this study from wells G8, G19 and G9 featured by Ca•Na-Cl water type with relative high TDS ranges (0.8–1.4, 1.3–1.6 and 1.5–2.8 g L⁻¹, respectively) and strontium concentrations (1.1–1.9, 4.9–7.1 and 7.3–11.6 mg L⁻¹, respectively), showing similarity with the geothermal system. Low temperature groundwater sampled in this study had Sr ∕ Cl mass ratios ranging from 2.4 × 10⁻⁴ to 1.6 × 10⁻², with higher ratios in deep groundwater (range: 9.4 × 10⁻⁴ to 1.3 × 10⁻², median: 3.7 × 10⁻³) compared to shallow groundwater (median: 3.1 × 10⁻³), and groundwater generally higher than saline surface water (range: 3.7 × 10⁻⁴ to 5.8 × 10⁻⁴, median: 3.9 × 10⁻⁴; Table S1).

Nitrate concentrations in groundwater range from 2.0 to 178.5 mg L⁻¹ (mean 90.1 mg L⁻¹) for shallow groundwater, and 2.0–952.1 mg L⁻¹ (mean 232.1 mg L⁻¹) for the deep groundwater, with most samples exceeding the WHO drinking water standard (50 mg L⁻¹).

5 Discussion

5.1 Groundwater isotopes and hydrochemistry as indicators of mixing processes

The Quaternary groundwater system in the Yang–Dai coastal river plain may be recharged by precipitation, irrigation return flow, river infiltration and lateral subsurface runoff (e.g., from mountain-front regions). Groundwater geochemical characteristics are then controlled by hydrogeological conditions and mixing processes, including mixing induced by extensive groundwater pumping, as well as natural mixing and water–rock interaction. It is evident from the geochemistry that mixing has occurred between groundwater and seawater in the coastal areas, as well as between normal/low temperature groundwater and geothermal water in the inland areas (e.g., near the Danihe geothermal field). Different sources of water are generally characterized by somewhat distinctive stable isotopic and hydrochemical compositions, allowing mixing calculations to aid understanding of the groundwater salinization and mixing processes, as discussed below.

Stable isotopes of O and H in groundwater and surface water fall on a best-fit regression line (dashed line in Fig. 3) with slope of δ²H = 4.4 × δ¹⁸O − 21.7, significantly lower than either the local or global meteoric water lines. Three processes are likely to be responsible for the measured range of isotopic compositions: (1) mixing between saline surface water (e.g., seawater or saline river water affected by tidal ingress) and fresher, meteoric-derived groundwater or surface water; (2) mixing between fresh meteoric-derived groundwater and saline thermal water; (3) evaporative enrichment of surface water and/or irrigation return flow, which may infiltrate groundwater in some areas. A subgroup of surface water samples (e.g., S1 to S3, S7 and S12; termed “brackish surface water”) show marine-like stable isotopic compositions and major ion compositions (Figs. 3 and 5). The “fresh” surface water samples (e.g., EC values < 1500 µS cm⁻¹) exhibit meteoric-like stable isotope compositions, with some samples (such as S9 and S10) showing clear evidence of evaporative enrichment in the form of higher δ²H and particularly, δ¹⁸O values (Fig. 3).

Figure 5 Relationship between chloride content and isotopic signature of different water samples as a means to differentiate mixing processes in the area. The data of thermal water are from Zeng (1991).

Fresh groundwater has depleted δ¹⁸O and δ²H values relative to seawater and show a clear meteoric origin, albeit with modification due to mixing. Theoretically, the mixing of meteoric-derived fresh groundwater and marine water should result in a straight mixing line connecting the two end-members; however, this is also complicated in the study area by the possible mixing with geothermal water. The thermal groundwater has distinctive stable isotopic and major ion composition (Han, 1988; Zeng, 1991), allowing these mixing processes to be partly delineated. Stable isotopes of thermal groundwater are more depleted than low-temperature groundwater (e.g., δ¹⁸O values of approximately −8 ‰, Fig. 8), indicating this likely originates from the mountainous areas to the north; Zeng (1991) estimated the elevation of the recharge area for the geothermal field to be from 1200 to 1500 m a.s.l. Based on a bivariate plot of δ¹⁸O vs. Cl⁻ with mixing lines and defined fresh and saline end-members, Fig. 5 shows the estimated degree of mixing between fresh groundwater, including shallow (G4) and deep (G25) groundwater end-members, and saline water, including seawater and geothermal end-members.

The two fresh end-members were selected to represent a range of different groundwater compositions and/or recharge sources, from shallow water that is impacted by infiltration of partially evaporated recharge (fresh but with enriched δ¹⁸O) to deeper groundwater unaffected by such enrichment (fresh and with relatively depleted δ¹⁸O). The narrower range and relatively enriched stable isotopes in shallow groundwater samples collected during the dry season compared with the wet season indicate some influence of seasonal recharge by either rainfall (fresh, with relatively depleted stable isotopes) or irrigation water subject to evaporative enrichment (more saline, with enriched stable isotopes and high nitrate concentrations; Currell et al., 2010) and/or surface water leakage. While there is overlap in the isotopic and hydrochemical compositions of shallow and deep groundwater (Figs. 3 and 4), this effect appears to only affect the shallow aquifer.

Based on Fig. 5, the shallow groundwater samples (e.g., G15, G10, G11, G14) collected from or around the Zaoyuan well field appear to be characterized by mixing between fresh meteoric water and seawater (plotting in the upper part of Fig. 5); while some deeper groundwater samples (e.g., G13, G2, G16, G14) collected from the coastal zone also appear to indicate mixing with seawater. Groundwater sampled relatively close to the geothermal field (e.g., G9, G19) shows compositions consistent with mixing between low-temperature fresh water and saline thermal water (lower part of Fig. 5). This is more evident in deep groundwater than shallow groundwater, which is consistent with mixing from below, as expected for the deep-source geothermal water. Other samples impacted by salinization show more ambiguous compositions between the various mixing lines, which may arise due to mixing with either seawater, geothermal water or a combination of both (e.g., G29).

Figure 6 Cross sections showing results obtained from application of geophysical resistivity method (employed in May 2004, data and methods described in Zuo, 2006).

The estimated mixing fraction (f_sw) of marine water for the shallow brackish groundwater ranges from 1.2 to 13.0 % and 2.6 to 6.0 % for the deep brackish groundwater. The highest fraction of 13 % was recorded in G10, located in the northern part of the Zaoyuan well field, which is located near a tidally impacted tributary of the Yang River (Fig. 1). Relatively higher fractions of marine water in relatively shallow samples (including those from the well field) compared to deeper samples may indicate a more “top down” salinization process, related to leakage of saline surface water through the riverbed, rather than “classic” lateral sea water intrusion, which typically causes salinization at deeper levels due to migration of a saltwater “wedge” (e.g., Werner et al., 2013); this is consistent with results of resistivity surveys conducted in the region (Fig. 6). The profile of chloride concentrations vs. depth indicates that salinization affects shallow and deep samples alike, with the most saline samples being relatively shallow wells in the Zaoyuan well field (Fig. 7).

Figure 7 Chloride concentration vs. well depth for groundwater samples.

Figure 8 Molar ratios of major ions vs. chloride concentrations for different water samples from the study area.

In general, brackish and fresh groundwater samples show distinctive major ion compositions, with the more saline water typically showing higher proportions of Na and Cl (Fig. 4). This contrasts with historic data collected from the Zaoyuan well field, which showed Ca-HCO₃ type water with Cl concentrations ranging from 130 to 170 mg L⁻¹. This provides additional evidence that the salinization in this area is largely due to marine water mixing. More Ca-dominated compositions are evident in the region near the geothermal well field further in-land (e.g., G5, G8, G19, G29 and G24), consistent with a component of salinization that is unrelated to marine water intrusion. Plots of ionic ratios of Na ∕ Cl and Mg ∕ Ca vs. Cl also reveal a subset of relatively saline deep groundwater samples which appear to evolve towards the geothermal-type signatures with increasing salinity (Fig. 8).

Stronger evidence of mixing of the geothermal water in the Quaternary aquifers (particularly deep groundwater) is provided by examining strontium concentrations in conjunction with chloride (Fig. 9). The geothermal water from Danihe geothermal field has much higher Sr concentrations (up to 89.8 mg L⁻¹) than seawater (5.4–6.5 mg L⁻¹ in this study), due to Sr-bearing minerals (i.e., celestite, strontianite) with Sr contents of 300–2000 mg kg⁻¹ present in the bedrock (Hebei Geology Survey, 1987). Groundwater sampled from near the geothermal field in this study has the highest Sr concentrations, e.g., G9 with Sr concentrations ranging from 7.4 to 11.6 mg L⁻¹ and G19 from 4.9 to 7.1 mg L⁻¹.

Figure 9 Chloride vs. strontium concentrations of different water samples. ZWF: Zaoyuan well field.

The plot of chloride vs. strontium concentrations (Fig. 9) shows that these samples and others (e.g., G16, G20, G27, G29) plot close to a mixing line between fresh low-temperature and saline thermal groundwater. Mass ratios of Sr ∕ Cl in these samples are also elevated relative to seawater by an order of magnitude or more (e.g., Sr ∕ Cl > 5.0 × 10⁻³, compared to 3.9 × 10⁻⁴ in seawater, Table S1). Other samples from closer to the coast (e.g., G4) also approach the thermal-to-low temperature mixing line, indicating probable input of thermal water. Samples collected from the Zaoyuan well field generally plot closer to the Sr ∕ Cl seawater mixing line (consistent with salinization largely due to marine water – Fig. 9); however, samples mostly plot slightly above the mixing line with additional Sr, which may indicate more widespread (but volumetrically minor) mixing with the thermal water in addition to seawater.

5.2 Anthropogenic pollution of groundwater

The occurrence of high nitrate (and possibly also sulfate) concentrations in groundwater in both coastal and in-land areas also indicates that anthropogenic pollution is an important process impacting groundwater quality and salinity (Fig. 10, Table 1). Seawater from Bohai Sea is heavily affected by nutrient contamination, showing ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations of 810 mg L⁻¹ in this study, and up to 1092 mg L⁻¹ in seawater further north of the bay near Dalian (Han et al., 2015), primarily due to wastewater discharge into the sea. The historic sampled ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentration of groundwater in the well field increased from 5.4 in May 1985 to 146.8–339.4 mg L⁻¹ in August 2010, while the concentration in seawater changed from 57.4 in May 1985 to 810.1 mg L⁻¹ in August 2010. A bivariate plot of Cl⁻ vs. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations in groundwater (Fig. 10) can thus be used to identify nitrate sources and mixing trends, including infiltration with contaminated seawater, and other on-land anthropogenic NO₃ sources (e.g., domestic or industrial wastewater discharge and/or ${\mathrm{NO}}_{\mathrm{3}}^{-}$-bearing fertilizer input through irrigation return flow).

Figure 10 Plot of chloride vs. nitrate concentrations in groundwater, with seawater and anthropogenic pollution mixing trajectories.

From this plot (Fig. 10) it appears that the major source of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in groundwater is on-land anthropogenic inputs rather than mixing with seawater, which would result in relatively large increases in Cl along with NO₃. Samples G10 and G15 (from the well field) are exceptions to this trend, showing clear mixing with nitrate-contaminated seawater. Deep groundwater (e.g., G9, G14) is also extensively contaminated with high ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations; this is likely associated with leakage from the surface via poorly constructed or abandoned wells – a problem of growing significance in China (see Han et al., 2016b; Currell and Han, 2017). According to one investigation by Zang et al. (2010), 14 of 21 pumping wells in the Zaoyuan well field have been abandoned due to poor water quality, and 307 pumping irrigation wells (occupied two-thirds of total pumping wells for irrigation) in the region have also been abandoned. Local authorities have, however, not implemented measures to deal with abandoned wells, meaning they are a future legacy contamination risk – e.g., by allowing surface runoff impacted by nitrate contamination to infiltrate down well annuli.

5.3 Hydrochemical evolution during salinization

A hydrogeochemical facies evolution diagram (HFE-D) proposed by Giménez-Forcada (2010) was used to analyze the geochemical evolution of groundwater during seawater intrusion and/or freshening phases (Fig. 11). In the coastal zone, the river water shows an obvious mixing trend between fresh and saline end-members. Some shallow groundwaters (e.g., G2, G4, G10, G13, G15) are also close to the mixing line between the surface-water end-members on this figure, indicating mixing with seawater without significant additional modification by typical water–rock interaction processes (e.g., ion exchange). Most brackish groundwaters (e.g., G11, G16, G17, G20, G25, G28, G29) have evolved in the series Ca-HCO₃ → Ca-Cl → Na-Cl, according to classic seawater intrusion. A relative depletion in Na (shown in lower than marine Na ∕ Cl ratios) and enrichment in Ca (shown as enriched Ca ∕ SO₄ ratios) is evident in groundwater with intermediate salinities (e.g., Figs. 8 and 11), indicating classic base exchange between Na and Ca during salinization (Appelo and Postma, 2005). Locally, certain brackish water samples (e.g., G1, G12, G26) appear to plot in the “freshening” part of the HFE diagram (potentially indicating slowing or reversal of salinization due to reduce in groundwater use), although these do not follow a conclusive trajectory. Water samples from the geothermal field (G5, G8, G9, and G19) plot in a particular corner of the HFE diagram away from other samples (being particularly Ca-rich); a result of their distinctive geochemical evolution during deep transport through the basement rocks at high temperatures.

Figure 11 Hydrogeochemical facies evolution diagram (HFE-D) for the collected water samples. Groundwater: gw.

5.4 Conceptual model of salinization and management implications

Coastal zones encompass the complex interaction among different water bodies (i.e., river water, seawater and groundwater). The interactions between surface water and groundwater in the Yang–Dai coastal river plain have generally been ignored in previous studies. However, the surface water chemistry data show that the distribution of saltwater has historically reached more than 10 km inland along the estuary of the Yang River, and approximately 4 km inland in the Dai River (Han, 1988). The relatively higher proportion of seawater-intrusion-derived salinity in shallow samples in this study, along with the evidence from resistivity surveys (Fig. 6, Zuo, 2006) indicate that intrusion by vertical leakage from these estuaries is therefore an important process. The hazard associated with this pathway in recent times has been reduced by the construction of a tidal dam, which now restricts seawater ingress along the Yang estuary to within 4 km of the coastline. This may alleviate salinization to an extent in future in the shallow aquifer by removing one of the salinization pathways; however, as described, there are multiple other salinization processes impacting the groundwater in the Quaternary aquifers of the region.

A conceptual model of the groundwater flow system in the Yang–Dai coastal river plain is summarized in Fig. 12. This model presents an advance on the previous understanding of the study area, by delineating four major processes responsible for groundwater salinization in this area. These are as follows: (1) seawater intrusion by lateral subsurface flow; (2) interaction between saline surface water and groundwater (e.g., vertical leakage of saline water from the river estuaries); (3) mixing between low-temperature groundwater and deep geothermal water; and, (4) irrigation return flow and associated anthropogenic contamination. Both the lateral and vertical intrusion of saline water are driven by the long-term overpumping of groundwater from fresh aquifers in the region. The irrigation return flow from local agriculture results from overirrigation of crops and is responsible for extensive nitrate pollution (up to 340 mg L⁻¹ ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in groundwater of this area), probably due to dissolution of fertilizers during infiltration. The somewhat enriched stable isotopes in shallow groundwater (more pronounced in the dry season) also indicate that such return flow may recharge water impacted by evaporative salinization into the aquifer. The geothermal water, with distinctive chemical composition (e.g., depleted stable isotopes, high TDS, Ca and Sr concentrations), is also demonstrated in this study to be a significant contributor to groundwater salinization, via upward mixing. The study area is therefore in a situation of unusual vulnerability, in the sense that it faces salinization threats simultaneously from lateral, downward and upward migration of saline water bodies.

Figure 12 Conceptual model of groundwater flow system in the Yang–Dai coastal river plain. Explanation: (1) land surface; (2) bedrock; (3) boundary between quaternary sediments and bedrock; (4) fault; (5) permeable fracture zone; (6) concentrated groundwater pumping zone; (7) pumping wells; (8) zone affected by upflow of geothermal fluids; (9) water table; (10) potential interface between fresh- and saltwater; (11) groundwater flow direction in Quaternary aquifers; (12) groundwater flow in bedrocks; (13) geothermal groundwater flow direction; (14) irrigation return flow.

According to drinking water standards and guidelines from the China Environmental Protection Authority (GB 5749-2006) and/or US EPA and WHO, chloride concentration in drinking water should not exceed 250 mg L⁻¹. At the salinity levels observed in this study – many samples impacted by salinization contain > 500 mg L⁻¹ of chloride (Table 1) – a large amount of groundwater is now or will soon be unsuitable for domestic usage, as well as irrigation or industrial utilization. So far, this has enhanced the scarcity of fresh water resources in this region, leading to a cycle of groundwater table decline → seawater intrusion → loss of available freshwater → increased pumping of remaining fresh water. If this cycle continues, it is likely to further degrade groundwater quality and restrict its usage in the future. Such a situation is typical of the coastal water resources “squeeze” highlighted by Michael et al. (2017). Alternative management strategies, such as restricting water usage in particular high-use sectors, such as agriculture, industry or tourism, that are based on a comprehensive assessment of the social, economic and environmental benefits and costs of these activities, warrants urgent and careful consideration.

6 Conclusions

Groundwater in the Quaternary aquifers of the Yang–Dai coastal river plain is an important water resource for agricultural irrigation, domestic use (including for tourism) and industrial activity. Extensive groundwater utilization has made the problem of groundwater salinization in this area increasingly prominent, resulting in the closure of wells in the area. Based on the analysis of hydrochemical and stable isotopic compositions of different water bodies, we delineated the key groundwater salinization processes. Seawater intrusion is the main process responsible for salinization in the coastal zone; however, this likely includes vertical saltwater infiltration along the riverbed into aquifers as well as lateral seawater intrusion caused by pumping for fresh groundwater at the Zaoyuan well field. The upward mixing of high TDS-geothermal water into the Quaternary aquifers is also evident, particularly through the use of stable isotope, chloride and strontium end-member mixing analysis. Additionally, significant nitrate pollution from the anthropogenic activities (e.g., agricultural irrigation return flow with dissolution of fertilizers) and local intrusion of heavily polluted seawater are also evident.

Groundwater salinization has become a prominent water environment problem in the coastal areas of northern China (Han et al., 2014, 2015, 2016a) and threatens to create further paucity of fresh water resources, which may prove a significant impediment to further social and economic development in these regions. Since the 1990s, the local government has begun to pay attention to the problem of seawater intrusion, and irrational exploitation of groundwater has been restricted in some cases. The Zaoyuan well field ceased to pump groundwater in 2007, while an anti-tide dam (designed to protect against tidal surge events) established in the Yang River estuary may also reduce saline intrusion in future. However, due to the significant lag time associated with groundwater systems, a response in terms of water quality may take time to emerge, and in the meantime the other salinization and pollution impacts documented here may continue to threaten water quality. In this regard, we recommend continued monitoring of groundwater quality and levels, and active programs to reduce input of anthropogenic contaminants such as nitrate from fertilizers, as well as appropriate well construction and decommissioning protocols to prevent contamination through preferential pathways.