The hydrogeochemical analysis and numerical simulation were conducted to explore the impact of land reclamation on submarine groundwater discharge and groundwater evolution in the seawater intrusion-retreat processes in eastern Jiangsu, China. The hydrogeochemical results indicated that the origin of salt in most groundwater is saltwater mixing, and weathering-dissolution of minerals enriches the salt in fresh groundwater. The evolution path of coastal groundwater in Jiangsu is Ca-Na-HCO3-Cl → Ca-Na-Cl-HCO3 → Na-Ca-Cl-HCO3 → Na-Mg-HCO3 → Na-Mg-HCO3-Cl → Na-Cl-HCO3 → Na-Cl, and the difference in groundwater evolution in different regions is obvious. The numerical simulation revealed that the land reclamation has blocked groundwater exchange and restrained the seawater intrusion. The sensitive analysis indicated that the permeability and the slope of the land reclamation have dominated the salinity distribution. And the aquifer heterogeneity promotes groundwater evolution and leads to the irregular salinity distribution of the groundwater. Moreover, the groundwater evolution in the shallow aquifer with the land reclamation shows that the residual seawater is gradually desalinating and discharged into the sea. The land reclamation is salinized by saltwater discharge and seawater intrusion together and then subsides into a stable saltwater wedge. The groundwater salinity distribution in Jiangsu is not fit seawater intrusion, which is caused by the coastline changes due to anthropogenic land reclamation and natural activities.

Seawater intrusion (SWI) is a major environmental issue in coastal areas with intensive anthropogenic activity. It deteriorates groundwater quality and decreases water availability, affecting industrial-agricultural production and resident living [14]. In recent years, the increasing groundwater exploitation due to population growth and urbanization in coastal areas has exacerbated SWI [58]. Therefore, investigating the factors and mechanisms related to SWI and groundwater evolution in intruded aquifers became important to prevent SWI and reduce its negative effects.

SWI and seawater retreat (SWR) in the coastal aquifers occur alternately [911], affecting the intrusion distance. However, the intruded seawater is restricted and accumulated in the layered heterogeneous aquifer. Additionally, the width of the freshwater-saltwater transition zone is larger compared with the homogeneous aquifer [12].

In addition to the heterogeneity of aquifers caused by natural faults and sediments, disturbances caused by anthropogenic activities may be more intense [13], and groundwater flow in the medium would be affected [1416]. Land reclamation (LR) is the anthropogenic construction aiming to expand the land area by extending the coast to the sea, which is widely found in coastal areas around the world, such as Holland, Japan, Saudi Arabia, Singapore, Korea, and China. LR increases the heterogeneity of the aquifer, the risk of natural disasters, and the ecological vulnerability in the coastal area [17]. It also reduces the carrying capacity of resources and intensifies the conflict between the economic-social system and the environment. Specifically, LR changes the hydrodynamic conditions of the ocean, causing a series of environmental issues, such as sedimentation, weakening of the wave-eliminating capacity in shallow offshore waters, and destructive effects of storm surges. LR has caused the disappearance of beaches, tidal flats, and natural wetlands, which has also destroyed coastal tourism resources and animal feeding grounds, leading to the endangerment of rare species. LR breaks the balance of interdependence between the sea and the land by creating a straight artificial coastline that reduces the area of the bay and estuary, blocking the river channel into the sea, increasing the magnitude and occurrence of floods. Overall, LR changes the characteristics of the water cycle between the surface water, groundwater, and seawater.

The methods to study SWI include numerical models [1820], analytical methods [2123], field hydrologic measurements [2426], microbial community researches [27, 28], geophysical methods [2931], and laboratory experiments [12, 32, 33]. Numerical simulation is an effective method for studying seawater intrusion. It can reflect the instantaneous and continuous state of groundwater flow and salt distribution in heterogeneous aquifers [32, 34, 35], mixing zone widths and seawater circulation rates.

The main objective of this paper is to explore the groundwater evolution and hydrochemical characteristics in intruded coastal aquifer with LR in eastern Jiangsu, China. Specifically, the study (a) builds numerical models to reveal the salinity distribution and groundwater flow in the seawater intrusion-retreat process, (b) explores the impact of the layered heterogeneity on groundwater evolution, and (c) performs a sensitivity analysis on LR, investigating the influence of factors, such as the LR permeability, thickness, and its sloping interface with the aquifer.

The Jiangsu Province is located on the west bank of the Yellow Sea, near the Yangtze River Delta. It belongs to the East Asian monsoon climate transition zone, with both continental and oceanic climate characteristics. The annual average temperature ranges from 13°C to 16°C. And the annual average precipitation is between 900 mm and 1100 mm, which concentrates from June to September and gradually increases from south to north.

The study area is located in the eastern coastal area of the Jiangsu Province, which is bounded between 32°00−35°5 N and 119°00−121°40 E. The regional aquifers can be divided into the phreatic aquifer group and the I−V confined aquifer group. The phreatic aquifer and the confined aquifer I are shallow pore water aquifer systems with no continuous aquiclude between them, which is dominated by brackish water or saline water. Among them, the phreatic aquifer group is composed of Holocene marine or sea-land alternate fine-grained sediments. The depth of the phreatic water level is generally 1 m−3 m, with the aquifer thickness range of 20 m−40 m. The desalination layer has an approximate thickness of 3 m and is located at the upper of the phreatic aquifer, containing brackish water. The confined aquifer group I is composed of loose sediments intersected by the late Pleistocene marine-terrestrial alternating facies. The depth of the confined aquifer I ranges from 40 m to 60 m, containing mainly brackish water. Shallow groundwater is mainly replenished by atmospheric precipitation and is closely linked to the rivers. And the leaching-infiltration of salty soil makes the shallow groundwater salinize. In contrast, the middeep groundwater is mostly freshwater, generally found in a semiclosed or closed state. So far, the precipitation from the 1950s has not been able to seep down to 33 m underground. Therefore, it is difficult for the middeep aquifer to receive recharge from precipitation or surface water [3638].

From 1128 to 1855, the Yellow River captured the Huai River and entered the Yellow Sea. As a result, a large amount of sediment quickly silted up to form the coastal plain of northern Jiangsu. Similar to the modern LR process, the coastline continued to move toward the sea (see Figure 1) [3941].

3.1. Water Sampling and Analysis

In April 2019, 59 shallow groundwater samples and two estuary seawater samples were collected from the coastal area of Jiangsu. They were collected from the nine sections, from north to south: Lianyungang I, Lianyungang II, Xiangshui, Binhai, Sheyang, Yancheng, Dongtai, Haian, and Rudong (see Figure 1). The collected samples were encapsulated in polyethylene plastic bottles after being filtered through a 0.45 μm membrane and then stored at 4°C until analyzed.

The total dissolved solids (TDS) were analyzed with a portable multiparameter water quality analyzer (YSI ProPlus, American YSI). HCO3 was tested by the titration method with phenolphthalein and standard HCl solutions. Major anion concentrations (F, Cl, Br, and SO42−) were analyzed using the ions chromatography (ICS-3000, American DIONEX). The cation concentrations (K+, Ca2+, Na+, Mg2+, and Sr) were measured with an inductively coupled plasma-atomic emission spectrometer (ICP-MS, American Thermo Fisher). The samples were acidified to pH2 with HNO3 for the cation test. The charge balance error in all sample analyses was less than 8%. The stable isotopic compositions (18O and 2H) were analyzed with an LGR (Los Gatos Research) liquid water isotope laser mass spectrometry. The values of δ18O and δ2H were calculated based on the Vienna Standard Mean Ocean Water (V-SMOW). The analytical precision of the long-term standard measurement of δ18O and δ2H were ±0.2‰ and ±0.6‰, respectively.

All groundwater samples were classified into four types based on TDS: freshwater (<1 g/L), brackish water (1–3 g/L), saline water (3–50 g/L), and brine (>50 g/L).

3.2. Scenario Definitions

We defined a set of scenarios to explore the salinity distribution and the groundwater evolution in the SWI and SWR processes (Table 1). The 2D vertical cross-sectional scenarios represent coastal homogeneous and layered heterogeneous unconfined aquifer (Figure 2 for a conceptual model). The aquifer setup assumes uniform flow perpendicular to the coast. The model domain has an extent of 200 m length and 42 m depth, which is a simplification of the regional unconfined aquifer. A no-flow boundary was defined on the upper and lower parts of the numerical model. The left-side freshwater boundary and the right-side seawater boundary were both set as a constant hydraulic head boundary, and the fixed salt concentration of the freshwater boundary and the seawater boundary is set to 0 kg/m3 and 35 kg/m3, respectively. The hydraulic head of the freshwater side is always higher than that of the seawater side, and the hydraulic gradient is 0.0075. These fall within a realistic range.

3.3. Numerical Simulations

SEAWAT is one of the most popular numerical codes for simulating variable-density groundwater flow and solute transport. It couples MODFLOW and MT3DMS with density-dependent terms in the governing equations. The physical processes of the seawater-freshwater interaction (FSI) and salinity distribution can be studied by SEAWAT. SEAWAT has been widely used for numerical simulations of SWI and calculation of the submarine groundwater discharge [4248]. The domain of the model is discretized into rectangular grids with 200m×1m×42m, which is referenced to the geological condition in the study area. The hydraulic parameters (hydraulic conductivity anisotropy, porosity, and dispersivity) were modified until the FSI and the stable saltwater wedge in the model are fit [4951] (see Table 2).

4.1. Occurrence of Water and Salt

Groundwater salinity is closely related to the origin of salt and its accumulation process. In coastal aquifers, the salt sources of the groundwater include the weathering-dissolution of minerals, SWI, saltwater mixing, and anthropogenic activities. The Br/Cl ratio is an important indicator to determine the salt source of the groundwater, especially for distinguishing between terrestrial evaporation origin and marine origin [52]. If the weathering-dissolution of evaporite is the dominant source of Cl in groundwater, then the Br/Cl ratio decreases with the increase of Cl concentration. If the Cl in groundwater is enriched by evaporation or seawater mixing, then the Br/Cl ratio remains constant with the increase of Cl concentration. The enriched Br in residual seawater or brine is caused by the greater solubility of bromide compared to that of NaCl. For fresh groundwater, the low Cl concentration with a large range of the Br/Cl ratio is caused by anthropogenic activities. Therefore, the Br sources in fresh groundwater cannot be accurately determined [53].

The Br concentration in the fresh groundwater ranges from 0.1285 mg/L to 227.2 mg/L, with an average value of 73.39 mg/L (Figure 3(a)). Most of the brackish water and saline water samples are distributed between the seawater Br/Cl line and the mixing line. These phenomena indicate that the Br source in groundwater is not only SWI but also the mixing of residual saltwater mixing in the study area. The process I (the Br/Cl ratio increases with the constant Cl concentration in fresh groundwater) shows that there are other sources of Br in groundwater besides seawater and residual saltwater in the study area. Process II (the Br/Cl ratio remains constant with increasing Cl concentration) reflects the influence of evaporation on groundwater salinity, while process III (the Br/Cl ratio decreases with increasing Cl concentration) reflects that groundwater is affected by residual saltwater mixing (Figure 3(b)). Thus, the salt sources in groundwater are SWI, residual saltwater mixing, and the weathering-dissolution of surrounding rock minerals in the study area. The weathering-dissolution of minerals is the dominant origin of salt in fresh groundwater.

The Na+/Cl ratio is a hydrogeochemical parameter that characterizes the degree of Na+ enrichment in groundwater. Additionally, the ratio can trace the origin of groundwater and the factors influencing salinity enrichment [54]. The Na+/Cl ratio in fresh groundwater increases, while the Cl concentration remains constant (Figure 4), indicating that silicate weathering is the dominant source of salt in fresh groundwater. For the brackish water and saline water, the Na+/Cl ratio shows a slight decrease along the seawater Na+/Cl line with the increase of Cl concentration. These phenomena indicate that the brackish and saline waters were mainly affected by saltwater mixing and SWI, and the Na+/Cl ratio further increases due to evaporation.

CaF2 is considered the dominant source of F in groundwater [26]. The increase of (Na++K+)/Ca2+ and HCO3 in groundwater promotes the release and enrichment of F in the surrounding rocks [55]. Sr in nature mainly exists in the soil and seawater. And the concentration of F and Sr in groundwater increases with the increase of seawater mixing, while the Sr/Cl ratio gradually decreases (Figure 5). The weathering-dissolution of minerals containing Sr or F leads to the enrichment of Sr/F in groundwater, especially for fresh groundwater.

The stable isotopic compositions (δ2H and δ18O) in coastal groundwater are related to recharge, mixing, and evaporation [26]. The groundwater is scattered near to the global meteoric water line (GMWL) from Craig [56] and the local meteoric water line (LMWL) of Jiangsu (Figure 6(a)), indicating that precipitation is the main recharge of groundwater. The LMWL of Jiangsu (δ2H=8.49δ18O+17.71, R2=0.9669, N=58) was provided by the International Atomic Energy Agency (IAEA). Compared with the GMWL (δ2H=8δ18O+10), the high slope and intercept of the LMWL reflect the humidity and rainy climate of the study area. These results indicate that the ocean moisture carried by the monsoon is the dominant source of precipitation due to this area close to the ocean. In addition, the stable isotopic compositions of groundwater reflect that the values of δ2H and δ18O gradually enrich with the increasing salinity (Table 3), and the groundwater samples are closer to the mixing line compared with the meteoric water line. These results indicate that the δ2H and δ18O of groundwater gradually enrich with the mixing of saltwater affected by evaporation.

4.2. Salinity Distribution and Groundwater Flow

Groundwater samples can be grouped based on the groundwater salinity distribution. For Lianyungang I, Xiangshui, and Sheyang, the groundwater salinity changes slightly without changes in the hydrochemical type. For Binhai and Yancheng, groundwater salinity increases first and then decreases with the hydrochemical type changing (brackish water → saline water → brackish water from sea to land). For Dongtai, Haian, and Rudong, the groundwater salinity gradually decreases with the hydrochemical type changing (saline water → brackish water → freshwater from sea to land), which is the common SWI characteristic. Finally, for Lianyungang II, the groundwater salinity increases with the hydrochemical type changing (freshwater → brackish water → saline water from sea to land).

The diversity of the groundwater salinity distribution indicates that the heterogeneity of the coastal aquifer is enormous perpendicular to the coast, and the SWI and SWR occur repeatedly. For example, the increase of groundwater salinity from sea to land indicates that saltwater is blocked by the low-permeability LR near the shoreline and remains in the original aquifer away from the modern coastline. Additionally, the groundwater in LR is freshwater or brackish water, indicating that the occurrence of SWI and SWR is restrained by geological conditions.

4.3. Hydrochemical Evolution

The Piper diagram is usually used to indicate the dominant ion compositions in water. The distribution characteristics of groundwater can be used to analyze groundwater evolution (Figure 7).

The coastal groundwater in Jiangsu can be divided into four types. The first category includes Lianyungang I, Xiangshui, and Sheyang, where freshwater, brackish water, and brackish water are distributed in the whole region, respectively. The second category includes Binhai and Yancheng, where brackish water, saline water, and brackish water are distributed in sequence from sea to land. The third category includes Dongtai, Haian, and Rudong, where saline water, brackish water, and freshwater are distributed in sequence from sea to land. Finally, freshwater and saline water are distributed sequentially from sea to land on the Lianyungang II section.

The cations of freshwater in Lianyungang I are dominated by Ca2+, followed by Na+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type of freshwater is Ca-Na-HCO3-Cl and Na-Ca-HCO3-Cl. The hydrochemical evolution of groundwater is Ca-Na-HCO3-Cl → Na-Ca-HCO3-Cl. The cations of brackish water in Xiangshui are dominated by Na+, followed by Mg2+. The anions are dominated by Cl, followed by HCO3. The hydrochemical type of brackish water is Na-Mg-HCO3-Cl and Na-Cl-HCO3. The hydrochemical evolution of groundwater is Na-Cl-HCO3 → Na-Mg-HCO3-Cl. The cations of brackish water in Sheyang are dominated by Na+, followed by Mg2+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type of brackish water is Na-Cl-HCO3, Na-Mg-HCO3-Cl, and Na-Mg-HCO3. The hydrochemical evolution of groundwater is Na-Cl-HCO3 → Na-Mg-HCO3-Cl → Na-Mg-HCO3.

The cations of brackish water in Binhai are Na+, followed by Mg2+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type is Na-Mg-HCO3-Cl, Na-Mg-HCO3, and Na-Cl-HCO3. The hydrochemical type of saline water is Na-Cl. The hydrochemical evolution of groundwater is Na-Mg-HCO3-Cl → Na-Mg-HCO3 → Na-Cl → Na-Cl-HCO3. The cation of brackish water in Yancheng is dominated by Na+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type is Na-HCO3-Cl and Na-Cl-HCO3. The hydrochemical type of saline water is Na-Cl. The hydrochemical evolution of groundwater is Na-HCO3-Cl → Na-Cl → Na-Cl-HCO3.

The hydrochemical type of freshwater in Dongtai is Na-Mg-HCO3. The cations of brackish water are dominated by Na+, followed by Mg2+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type is Na-Mg-HCO3-Cl and Na-Cl-HCO3. The hydrochemical evolution of groundwater is Na-Cl-HCO3 → Na-Mg-HCO3-Cl → Na-Mg-HCO3. The cations of freshwater in Haian are dominated by Mg2+ and Na+. The anions are dominated by HCO3. The hydrochemical type is Mg-Na-HCO3. The cations of brackish are dominated by Na+, followed by Mg2+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type is Na-Mg-HCO3-Cl and Na-Cl-HCO3. The hydrochemical type of saline water is Na-Cl. The hydrochemical evolution of groundwater is Na-Cl → Na-Cl-HCO3 → Na-Mg-HCO3-Cl → Mg-Na-HCO3. The hydrochemical type of freshwater in Rudong is Ca-Mg-HCO3. The cations of brackish water are dominated by Na+, followed by Mg2+. The anions are dominated by HCO3, followed by Cl. The hydrochemical type is Na-Mg-HCO3-Cl and Na-Cl-HCO3. The hydrochemical type of saline water is Na-Cl. The hydrochemical evolution of groundwater is Na-Cl → Na-Cl-HCO3 → Na-Mg-HCO3-Cl → Ca-Mg-HCO3.

The cations of freshwater in Lianyungang II are dominated by Na+, followed by Ca2+. The anions are dominated by Cl, followed by HCO3. The hydrochemical type is Na-Ca-HCO3-Cl, Na-Ca-Cl-HCO3, and Na-Ca-Cl. The hydrochemical type of saline water is Na-Cl. The hydrochemical evolution of groundwater is Na-Ca-HCO3-Cl → Na-Ca-Cl-HCO3 → Na-Ca-Cl → Na-Cl.

For the coastal groundwater in Jiangsu, the hydrochemical types of freshwater are Ca-Na-HCO3-Cl, Ca-Na-Cl-HCO3, Na-Ca-Cl-HCO3, and Na-Mg-HCO3. The hydrochemical type of brackish water is Na-Mg-HCO3, Na-Mg-HCO3-Cl, and Na-Cl-HCO3. The hydrochemical type of saline water is Na-Cl. Thus, the complete path of groundwater evolution is Ca-Na-HCO3-Cl → Ca-Na-Cl-HCO3 → Na-Ca-Cl-HCO3 → Na-Mg-HCO3 → Na-Mg-HCO3-Cl → Na-Cl-HCO3 → Na-Cl.

4.4. Numerical Simulation of Groundwater Evolution in the Intruded Aquifer with LR

4.4.1. SWR in Simplified Homogeneous Aquifers of the Jiangsu Area

In the coastal area with LR where SWI has occurred, the intruded seawater will be sealed within the aquifer by low-permeability LR. Under the influence of submarine groundwater discharge (SGD), the residual seawater in the intruded aquifer is continuously diluted and retreated, finally forming a stable saltwater wedge. For further research, two numerical experiments of the intruded homogeneous unconfined aquifer (hydraulic conductivity, K=1m/d) were simulated to explore the rate and path of SWR and salt migration. The experiments were set without LR and with the uninvaded LR (K of 0.1 m/d, LR slope of 45°, and bottom thickness of 21 m). The hydraulic head of the freshwater side is higher than that of the seawater side (fixed hydraulic gradient of 0.0075).

The numerical experiment of SWR in the intruded homogeneous aquifer (K=1m/d) without LR, and the hydraulic gradient of 0.0075 was simulated as the control experiment (Figure 8(a)). The density of seawater is greater than that of freshwater. Therefore, freshwater is discharged from the upper of the aquifer, forming a saltwater wedge. When it reaches 28,000 days, the shape and location of the saltwater wedge formed by SWR are consistent with that of the saltwater wedge formed by SWI in the uninvaded aquifer with the same conditions.

In the numerical experiment of the intruded aquifer with the LR (Figure 8(b)), the SWR on the freshwater side, the saltwater discharge, and the SWI on the seawater side occur simultaneously. As the simulation runs, seawater in the intruded aquifer and from the seawater side is discharged into the LR from the bottom and fills the LR from bottom to top. In the intruded aquifer, residual seawater is continuously discharged from the bottom and the vertical saltwater wedge becomes thinner. When the simulation runs approximately 4,200 days, the FSI in the intruded aquifer and the LR are consistent. Then, the saltwater thickness in the aquifer continues to be become thin, while the saltwater thickness in the LR increases firstly and then decreases. Ultimately, the saltwater wedge recedes into a stable state after 60,000 days of simulation. In conclusion, the low-permeability LR hinders the SWR and changes the salinity distribution of shallow groundwater along with the groundwater flow.

4.4.2. Extended Research on SWR in Layered Heterogeneous Aquifers

The geological conditions have a significant influence on the SWR and salt migration in the aquifer. For further research, a series of numerical experiments of varying coastal aquifer structures have been simulated to explore the rate and path of SWR and salt migration. These simulations include (a) the intruded layered heterogeneous unconfined aquifer (KT=1m/d, KM=0.1m/d, and KB=1m/d) without LR and (b) the intruded layered heterogeneous aquifers (KT=1m/d, KM=0.1m/d, and KB=1m/d; KT=0.1m/d, KM=1m/d, and Kb=0.1m/d; KT=1m/d, KM=0.1m/d, and KB=0.01m/d; and KT=0.01m/d, KM=0.1m/d, and KB=1m/d) with the uninvaded LR (the LR slope is 45°, its thickness is 21 m, and permeability is 0.1 m/d). The fixed hydraulic gradient is 0.0075 (the hydraulic head of fresh groundwater is higher than that of the seawater).

In the numerical experiment of the intruded layered heterogeneous aquifer (KT=1m/d, KM=0.1m/d, and KB=1m/d) without LR (Figure 9(b)). The fixed hydraulic gradient is 0.0075 (the hydraulic head of fresh groundwater is higher). The top and bottom layers with high permeability (K=1m/d) take the lead in SWR. Afterward, the SWR in the top layer is significantly faster than that of the bottom layer. The SWR in the middle low-permeability layer (K=0.1m/d) has a hysteresis, and the saltwater wedge in the whole aquifer is no longer smooth. In addition, there is an obvious desalination area in the bottom layer due to its permeability and the density difference between saltwater and freshwater. Compared with the homogeneous aquifer, SWR in the layered heterogeneous aquifer is blocked, and saltwater is discharged into the seawater side through channels that promote groundwater flow. This will form diversified salt distribution characteristics.

In the simulation of the intruded layered heterogeneous aquifer (KT=1m/d, KM=0.1m/d, and KB=1m/d) with the uninvaded LR (K=0.1m/d) (Figure 10(a)), SWR on the freshwater side and SWI on the seawater side occur simultaneously. On the freshwater side, the SWR is similar to the above-mentioned layered aquifer without LR. On the seawater side, saltwater discharge in the intruded aquifer and SWR occur simultaneously and the LR is filled by saltwater from the bottom. The saltwater in the top layer near the LR gradually retreats inland due to the push of freshwater. Compared with the layered heterogeneous aquifer without LR, LR promotes SWR vertically while inhibiting the saltwater wedge migration in the bottom layer.

In the simulation of the intruded layered heterogeneous aquifer (KT=0.1m/d, KM=1m/d, and KB=0.1m/d) with the uninvaded LR (K=0.1m/d) (Figure 10(b)), SWR on the freshwater side and SWI on the seawater side occur simultaneously. On the seawater side, the saltwater discharge from the middle layer (K=1m/d) and the seawater intrusion occur simultaneously, and the LR is filled with saltwater from the bottom. The saltwater in the top layer near the LR gradually retreats inland due to fresh groundwater flow. On the freshwater side, SWR in the middle layer (K=1m/d) occurs firstly, followed by the top layer (K=0.1m/d). Compared with the layered heterogeneous aquifer (KT=1m/d, KM=0.1m/d, and KB=1m/d), the rate of seawater retreat is much slower, and fresh groundwater is mainly discharged from the middle layer with high permeability (K=1m/d).

In the simulation of the intruded layered heterogeneous aquifer (KT=1m/d, KM=0.1m/d, and KB=0.01m/d) with the uninvaded LR (K=0.1m/d) (Figure 10(c)), the saltwater migration on the seawater side is similar to the above-mentioned layered heterogeneous aquifer with LR. Compared with the homogeneous aquifers and the above-layered heterogeneous aquifers, SWR concentrates on the top layer, the saltwater recedes slowly, and the driving force of freshwater discharge weakens. Thus, the time that it takes to salinize the LR with saltwater discharge and seawater intrusion is longer.

In the simulation of the intruded layered heterogeneous aquifer (KT=0.01m/d, KM=0.1m/d, and KB=1m/d) with the uninvaded LR (Figure 10(d)), the saltwater migration on the seawater side is similar to the above-mentioned layered aquifer with LR. On the freshwater side, SWR in the bottom layer (K=1m/d) occurs firstly, followed by the middle layer (K=0.1m/d). In contrast, the saltwater in the middle layer (K=0.1m/d) is diluted and retreated firstly, followed by the bottom layer (K=1m/d). The saltwater in the top layer is diluted from the bottom due to the density difference.

It takes more time for the saltwater wedge to reach a stable state, and the seawater intrusion distance becomes smaller compared with the aquifer without LR (Figure 11). The smaller permeability of the aquifer, the more time it takes for the saltwater wedge to reach a stable state, especially for the aquifer with the low-permeability medium at the bottom. And the bottom permeability of the aquifer and the SGD has affected the intrusion distance; the smaller permeability and the greater SGD, the smaller the intrusion distance.

When the hydraulic gradient is fixed, the residual seawater in the intruded aquifer continuously discharges until the interface of saltwater and freshwater balances. And the distance of the saltwater wedge toe decreases firstly and then increases; finally, the saltwater wedge is consistent with that of SWI under the same structure of the layered uninvaded aquifer with LR.

4.4.3. Sensitivity Analysis on the Land Reclamation

In the LR activities, the slope of the interface between aquifer and LR is affected by the coastal topography, while the permeability and the thickness of the LR are dominated by filled materials and the range of LR activities, respectively. These will affect the rate and the path of SWR and the salt migration in the intruded aquifer. For further research, three sets of numerical experiments of the intruded layered heterogeneous aquifer of type I were simulated, with the uninvaded LR and the fixed hydraulic gradient of 0.0075 (the hydraulic head of fresh groundwater is higher). The experiments explored the influence of the LR slope (30°, 45°, and 60°), permeability (0.1 m/d, 0.01 m/d, and 0.001 m/d), and thickness (1 m, 21 m, and 42 m) on the rate and the path of seawater intrusion-retreat and salt migration.

The SWR on the freshwater side and the SWI on the seawater side occur simultaneously. In the simulation experiment of the layered heterogeneous aquifer with the LR (Figures 10(a) and 12), as the LR slope increases, SWR is blocked and the residual seawater will stay longer. The slope increases the refraction angle of saltwater passing through the interface between the aquifer and LR, which salinizes the shallow groundwater in the LR for a longer time.

As the permeability of the LR decreases and the bottom thickness increases, the rate of SWR gradually decreases (Figures 10(a), 13, and 14). Compared with the thickness of LR, permeability plays a more dominant role in blocking the saltwater discharge and SWI on the seawater side, and the salinization of the LR is also prevented. The thickness difference of the LR slightly changes the salt distribution in the seawater intrusion-retreat, while the permeability significantly affects the salt distribution and leads to the diversity of groundwater evolution along with the groundwater flow. In addition, the slope of LR also affects the salt distribution, hinders the groundwater hydrochemical evolution, and leads to the diversity of groundwater hydrochemical types along with the groundwater flow.

The smaller the slope of the LR, the more time it takes for the saltwater wedge to reach a stable state and the greater the intrusion distance (Figure 15). And the smaller the permeability of the LR, the more time it takes for the saltwater wedge to reach a stable state and the smaller the intrusion distance. While the smaller the thickness of the LR, the less time it takes for the saltwater wedge to reach a stable state, and the intrusion distance decreases firstly and then increases with the increasing thickness of LR.

In the groundwater system with a complex aquifer structure, the rate and degree of groundwater exchange dominate the groundwater evolution. When the groundwater flow and the saltwater discharge have been blocked, the rate and extent of groundwater evolution are larger. In contrast, the complexity of the aquifer structure in the groundwater system leads to an irregular distribution of groundwater hydrochemical types in the vertical or horizontal directions.

4.5. Insights of the Experiment on Coastal Groundwater Evolution and Its Application in Reality

In the numerical simulation experiments of the intruded homogeneous unconfined aquifer with the LR, residual seawater in the aquifer is continuously discharged into the LR from the bottom. The wedge freshwater discharge channel is formed in the upper due to the density difference. In addition, the vertical thickness of the freshwater channel gradually increases. While in the LR, the residual seawater discharge and the SWI have continuously filled the LR from bottom to top, the LR is salinized firstly and then desalinized. Ultimately, saltwater recedes into a stable wedge in the aquifer with the LR system. Compared with the aquifer without LR, the LR blocks the seawater retreat, changes salt migration, and promotes groundwater evolution. In this process, the groundwater type in the shallow aquifer from the sea to the land gradually evolves as freshwater → saline water → freshwater → saline water → brackish water/freshwater → brackish water → saline water/brackish water → freshwater → saline water → brackish water/freshwater → freshwater → all freshwater. Among the nine sections from the sea to the land in Jiangsu, the section Lianyungang II corresponds to the salinity distribution of freshwater → saline water. Sections Binghai and Yancheng correspond to the distribution of freshwater → saline water → brackish water/freshwater and brackish water → saline water/brackish water → freshwater. Sections Dongtai, Haian, and Rudong correspond to the salt distribution of saline water → brackish water/freshwater → freshwater. Sections Lianyungang I, Xiangshui, and Sheyang correspond to full freshwater.

In the numerical experiments of the intruded layered heterogeneous aquifer with LR, permeability is the dominant control of the residual seawater distribution, followed by the density difference between seawater and freshwater. The upper layer with high permeability promotes the residual seawater discharge in the aquifer. In contrast, for the aquifer with LR, the groundwater flow velocity is significantly affected. The salt distribution is dominated by the slope and the permeability of the LR, which leads to the diversity of groundwater evolution along with the groundwater flow.

Regional high-salinity groundwater could be caused by leaching-infiltration of the salty soil. In addition, the poor permeability of the regional aquifer can prevent the occurrence of SWI. These factors lead to the irregular distribution of groundwater salinity, and it agrees with the several salinity distributions mentioned previously.

Eastern Jiangsu has experienced several natural and anthropogenic coastline expansion events over the years. Compared with the inevitable natural change of the coastline, the anthropogenic LR has led to great-scale ecological damage, especially for the submarine groundwater discharge, the groundwater purification, and the groundwater evolution in the land-sea interaction zone.

The selected ions (Br, Na+, Sr, F, and Cl) combined with the stable isotopic compositions (18O and 2H) in groundwater indicate that precipitation is the main recharge source of groundwater and the salt in most groundwater is dominated by saltwater mixing, while the weathering-dissolution of minerals promotes the salt enrichment in fresh groundwater. The groundwater salinity distribution reveals that groundwater in the study area can be divided into four types, and the residual seawater is restricted in the aquifers away from modern coastlines; meanwhile, modern SWI is also restrained by the LR.

During SWR and SWI, the LR blocks the groundwater flow and reduces the intrusion distance. Among them, the permeability and the slope of the LR dominate the salinity distribution. For the aquifer with the LR, the residual seawater is discharged from the bottom, and the LR is filled by saltwater due to the residual seawater discharge and seawater intrusion and then desalinized, eventually receding into a stable saltwater wedge. For the layered heterogeneous aquifer, permeability is the dominant control of the residual seawater distribution, followed by the density difference between seawater and freshwater. The high-permeability upper layer promotes the residual seawater discharge in the aquifer. In groundwater systems, complex aquifer structure hinders groundwater exchange and promotes the rate and extent of groundwater evolution. Meanwhile, the complex aquifer structure in the groundwater system leads to the irregular distribution of groundwater hydrochemical types in the vertical or horizontal directions.

The simulation results indicate that groundwater flow and salt migration conform to certain rules. Groundwater evolution in the shallow aquifer is freshwater → saline water (from the sea to the land) → freshwater → saline water (some saline water changing into brackish water/freshwater) → brackish water → saline water (some saline water changing into brackish water/freshwater) → saline water → brackish water/freshwater → freshwater. The four type groups in the eastern coastal area of Jiangsu are consistent with the groundwater evolution in the aquifer with LR.

The data are listed in the Supplementary Materials.

The authors declare that they have no conflicts of interest.

This work was funded by the National Key Research and Development Program of China (Nos. 2019YFC1804304 and 2016YFC0402801) and the Natural Science Foundation of China (NSFC Grants 41706067 and 42072274).

Supplementary data