Abstract

Detailed vertical profiles of Cl in porewaters through the aquitard-aquifer system were used to yield solute transport mechanism and build a conceptual model regarding evolution processes and transport time of natural tracer migration in North Jiangsu coastal plain, China. One-dimensional vertical simulated models of Cl profiles illustrate that diffusion appeared to be the dominant solute transport mechanism in the aquitard-aquifer system. A downward groundwater flow did not improve the fitness between simulated and measured values. Several simulated models were constructed and suggested that the evolution of the Cl profiles is mainly ascribed to the introduction of seawater and freshwater of transgression-regression to the first confined aquifer and the upper boundary. Groundwater in the first confined aquifer recharged by the Late Pleistocene glacial meltwater (25–15 ka BP) was supported in response to the low Cl concentrations. The shallow groundwater in the first confined aquifer and porewater with high salt were attributable to the Holocene seawater intrusion. These timeframes were also consistent favorably with the results of previous studies into the palaeohydrology of the study area.

1. Introduction

Aquitards frequently seem to have the ability to contain high salinity relative to aquifers [1, 2]. Many researches have been conducted on the thicker, nonfractured clay-rich and shale units with the bulk hydraulic conductivity (K) below 1 × 10–10 m/s [36]. Under such conditions, solute transport in the aquitards is demonstrated to be dominated by molecular diffusion [7, 8].

Shaw and Hendry [9] suggest that the thickness of clay-rich aquitard was requested to be >60 m in order to avoid the occurrence of the initial interferences among the advective-diffusive solute profiles. However, the water-bearing units (aquifers or sand streaks) with various thicknesses were often interspersed within aquitards [6, 1012] and particularly marked in the late Quaternary sediments of coastal plains in China due to the transgression-regression effect [13]. Critically, these occurrences could cause perturbations in transport paths and partial advection solute transport in porewater profiles and lead to problems with the interpretation of the palaeohydrogeological information [12]. Despite the potential for significant advective migration during tracer profiles, diffusion was determined as the pure transport mechanism controlling solute transport in the aquitard-aquifer system based on the researches of Kuang et al. [14, 15]. Geochemical tracers were used to define and constrain long-term transport mechanisms and preserve a historical record of the major palaeohydrologic events at the aquitard scale [6, 1012, 16]. Although solute transport in the aquitards had been conducted in many studies, the porewater transport mechanism and processes in clay-rich aquitards interspersed with the cooccurrence of multiple aquifers of coastal area need more studies to support and further confirm.

North Jiangsu coastal plain (NJCP) is located in the eastern part of China, adjacent to the south Yellow Sea (Figure 1). This study area belongs to the part of the new Silk Road, where groundwater resource plays a pivotal role in promoting social and economic development. However, groundwater in this region has been mostly seriously affected by seawater intrusion because of the Quaternary transgressions [17, 18]. To ensure the sustainable development of the new Silk Road, the management of the groundwater resource should be given considerable attention. Most of the previous studies in this area have focused on the detailed investigation of the change and evolution of regional groundwater quality [19, 20]. Additionally, saline groundwater in the aquitards and aquifers has been studied to investigate the hydrogeochemical characteristics and chemical evolution processes [21]. Nevertheless, an understanding of the saline transport mechanism and processes in the aquitard-aquifer system is still very limited in such coastal areas.

The aim of this work is to investigate the solute transport mechanism and evolution processes in the aquitard-aquifer system and provide additional palaeohydrogeologic information using Cl concentration vertical profiles. The specific objectives of this study are to (1) obtain the effective diffusion coefficients () and of the studied vertical deposition profiles on the basis of laboratory experiments, (2) apply 1D vertical transport model of the Cl profiles to gain insight into the dominant solute transport mechanism in the aquitard-aquifer system, and (3) explore a more detailed synthesis of the time of paleohydrologic and formative geologic events for tracer profiles in the aquitard-aquifer system.

2. General Description of the Field Site

2.1. Hydrogeological Setting

The current study is conducted on the Late Pleistocene and Holocene clay-rich deposits in NJCP, a part of China’s Eastern Plain (Figure 1). The NJCP has a continental and maritime climate and the average annual temperature of 13–16°C. Annual precipitation averaged approximately 800–1,200 mm and nearly 30–60% of the annual precipitation falls between June and September [25]. The annual evaporation is about 900–1,050 mm.

The depositional facies of the deposits change from west to east, progressing from alluvium to proluvial sediments and then marine sediments on the coastal plain, mainly because of the development of the Huai River, the Yellow River, and the sea level change [23, 26, 27]. A gradual and vertical Quaternary sediment shift is observable from the single continental alluvial to the transitional sediments between land and marine facies [17, 18]. The depositional facies of the research area can be classified into three groups: continental facies, marine facies, and transitional facies. Quaternary continental facies mainly refer to fluvial alluvial facies and fluvial-lacustrine alluvial facies [17, 18]. Marine facies are composed of littoral facies and shallow sea facies. Transitional facies occurred along the coast consist of estuarine, lagoon, intertidal, and residual seawater zone.

The major aquifer is composed of the multilayered aquifer groups on the east coast (Figure 2). The phreatic aquifer and the first, second, and third confined aquifers are recognized in this study area and composed of Quaternary sediments with a thickness of 10–50 m [19]. These aquifers consist of sandy gravel, medium-fine sand, and fine sand and are separated by silt- and clay-dominated aquitards. The aquifers are characterized by a complex multilayered framework due to the special geographical position and complex climate conditions (Figure 2).

The burial depth of the phreatic aquifer and the first confined aquifer are generally less than 5 and 60 m, respectively. Groundwater in the phreatic aquifer and the first confined aquifer corresponding to the Holocene and the Late Pleistocene formations is mostly saline with total dissolved solids (TDS) generally over 3 g/L. Freshwater can be found in the bottom of the first confined aquifer. The TDSs of saline water exhibit an increasing trend from the west to east. Groundwater in these aquifers belongs to Cl-Na water type, with a small amount of HCO3-Cl-Ca·Na type, which was derived from the Holocene seawater [20] and then mixed with the recharge of precipitation and irrigation water. Evaporation and exploitation are the main discharge pathways. Groundwater exploitation is limited due to the high salinity.

The burial depths of the second and third confined aquifers were usually shallower than 140 and 250 m, respectively. Groundwater in the second confined aquifer illustrated the chemical signature of fresh water and was trapped in the fluvial deposits during the Middle Pleistocene, which was possibly recharged during a colder period in the Late Pleistocene [28]. The third confined aquifer groundwater is mainly characterized by less saline with TDS below 2.0 g/L, except some areas in the northeast with TDS values extremely higher than 2.0 g/L [19, 20]. The salinity of some groundwater in the third aquifer was attributed to the entrapped relict seawater of the Late Pleistocene. The hydrochemical compositions of groundwater in the second and third confined aquifer are mostly HCO3-Na, HCO3·Cl-Na·Ca, Cl-Na, and Cl-Ca·Na. Groundwater runoff conditions are poor and the water circulation is slow gradually, which are current targets for exploitation and major water resource for local residents. Groundwater ages in the second and third aquifers were assumed to be around 30 ka BP on the basis of carbon isotopic dating (14C) [20, 28]. As a result of the lack of clay aquitards locally, the leakage recharge from the upper shallow groundwater is possibly regarded as the main sources for deep confined groundwater during groundwater exploitation, and the corresponding discharge is artificial extraction.

2.2. Transgressions and Regressions

It is believed that this study area has undergone many times of transgressive-regressive processes since Quaternary [29, 30]. In the early Middle Pleistocene, climate warming, strong surface runoff, and enhancement of river erosion caused the deposition of the clayey sediment interbedded with silt. Subsequently, the first transgression occurred as a result of sea level rise, and the range and degree were smaller.

The scope of the two other transgressions in the Late Pleistocene (Marine Isotope Stage (MIS) 5, about 110–70 ka BP; MIS 3, about 40–25 ka BP) expanded gradually [31]. These two transgressions corresponded to the depth of approximately 30–90 m in the study region. The shoreline of the second transgression (MIS 5) reached the above 5–7 m of the modern coastline, and all of the study area suffered from marine environment [29]. The transgression in the MIS 3 was relatively small, and the deposition environment might be subaerial exposure [32]. The sea level in MIS 3 seemed to not reach the location of the modern sea level in the study area and generally agreed that the coastline was close to or below 10–20 m of the modern shoreline [29]. In other words, the sea level regressed from this region with the time of 70–10 ka BP (MIS 2–4), and the corresponding sediments probably suffered from weathering denudation [33]. In the Holocene (MIS 1), the sea level rose rapidly and large-scale seawater intrusion occurred, resulting in an interactive marine and terrestrial deposit with a whole transgression-to-regression succession [22, 34]. It was not until AD 1,128 that the coastline then moved eastward rapidly and seawater gradually retreated from this region [35].

2.3. Boreholes Description

Two stratigraphic boreholes were drilled about 10 km apart (Figure 1). The first borehole (SY1) had a depth of 250 m and was located approximately 10 km west of the southern Yellow Sea shore (33.86° N, 120.43° E). In this study, the upper section with the depth of 2–182 m was conducted. The second borehole had a depth of 120 m and was located approximately 20 km west of the southern Yellow Sea shore (33.80° N, 120.32° E). The stratigraphic columns of the two boreholes were presented in Figure 3 reference to the field logs. 1-2 m below ground was considered to be the oxidized and fractured layer due to the fluctuation of water table, and the underlying deposits were regarded as unoxidized zone. There was a lack of downward water flow in the oxidized layer because of its fractured features and discontinuity [36]. The oxidized layer was outside the range of this study. Although the two boreholes were very close to each other, the marine layers have different degrees of inclination toward the sea because of neotectonic activity [13].

Three transgressive events were confirmed to have occurred in the study area and were named, from the earliest to the latest: the Asterorotalia transgression, the False Rotifer transgression, and the Winding Worms transgression. Species such as Ammonia annectens, A. beccarii var., and Rosalina bradyi were found in the surficial aquitard (grey-yellow) with a depth of 0–19.9 m in borehole SY1 and 0–18.9 m in borehole SY2, suggesting that these depths belonged to MIS 1 period. Correspondingly, the occurrence of a large number of Pseudorotalia schroeteriana indicates the underlying aquitard at depth of 23.8–33.6 m in SY1 and 23.3–54.2 m in SY2 formed in MIS 3. Deposits at depth of 41.5–86.5 m in SY1 and 61.1–85.6 m in SY2 developed during the MIS 5 period and were attributable to the presence of the warm species Pseudorotalia indopacifica and Asterorotalia pulchella.

3. Material and Methods

3.1. Borehole Sampling and Analysis

In this study, core samples were collected from 2 to 182 m in borehole SY1 between January and February 2013 using a rotary drill. Additional core samples were collected from 10.6 m to 120 m in borehole SY2 in December 2014. For these study profiles, subsamples were taken every 2–5 m at the SY1 and SY2 drilling sites. To ensure that no drilling fluid contaminated the collected samples, immediately after retrieval and before packing, the outer 2-3 cm of the samples were removed and discarded. Sample was packaged in a sealable aluminum barrel (in length of 20 cm and diameter of 8 cm) and sealed with wax, which was prepared for chemical analysis. Samples were collected on-site by 10 mL glass bottles to extract porewater for stable isotopic analysis and were sealed using raw adhesive tape immediately after sampling to prevent fractionation caused by evaporation. The packaged borehole samples were placed in coolers (at approximately 4-5°C) with ice bags in the field for transport to the School of Environmental Studies of the China University of Geoscience, where samples were stored at 4-5°C prior to analysis, to minimize the growth of microorganisms.

The high-pressure mechanical squeezer (HPMS) using nitrogen was employed on the sealed core samples for porewater extraction. The HPMS was designed and developed based on the ex situ squeezing device described in detail by Li et al. [37]. The application of a pressure-controlled noble gas to a sample forced porewater to separate from the sediment. The operational pressure of the HPMS was controlled within an allowable range of 0–8 MPa, with an analytical precision of 0.2 MPa. The outermost 1-2 cm of the prepared core samples was scraped to discard any materials that could have been altered by exposure to the atmosphere (oxidation). Approximately 1,000 g soil samples were conducted to ensure that enough porewater was provided for analysis. To avoid overconsolidation or destruction of the clay-pore system, the applied stress was increased to 8 MPa gradually rather than in a single step. First, a small stress of approximately 1 MPa was initially applied to expel most of the air in the sample chamber and to ensure that the sample was bedded. Second, the applied stress was increased by 1  MPa every 4-5 hours during the day. Finally, core samples were compressed and consolidated under pressure for 3-4 days. Meanwhile, porewater flowed out from the center hole of the bottom endplate through 0.45 μm filter paper and was collected in a clean plastic bottle. The water samples were weighed and filtered through a 0.45 μm membrane and immediately stored in a refrigerator at 4°C before being analyzed.

For stable isotopes analysis, porewater was extracted from the collected core samples using vacuum distillation, reducing the influence of air during the porewater collection process. The final distillation temperature keeps 120°C for 6-7 h. All extracted water samples were kept in a refrigerator at around 4°C.

Porewater samples of borehole SY1 () and SY2 () were analyzed for Cl and Br concentration analysis using ion chromatography (ICS-1100, Dionex, Sunnyvale, CA, USA) at the School of Environmental Studies, China University of Geosciences. The analytical precision for Cl and Br concentration was greater than 0.01 mg/L.

Stable isotope analysis in borehole SY1 and SY2 was conducted at the State Key Laboratory of Water Resources and Hydropower Engineering Science of Wuhan University and the Laboratory of Geological Survey Institute, China University of Geosciences in Wuhan, respectively. These isotope samples of SY1 and SY2 boreholes were analyzed with mass gas isotopic ratio mass spectrometry (MAT 253, Thermo Fisher Scientific, Waltham, MA, USA, with an analytical precision of ±0.2 and ±2, resp.) and a liquid water isotope analyzer (LGR, IWA–45EP, USA, with an analytical precision of ±0.1 and ±0.5, resp.), with values reported relative to Vienna Standard Mean Ocean Water (V–SMOW).

3.2. Laboratory Test
3.2.1. Porosity

The porosities of the soil samples from the saturated zone were determined. To avoid samples being in an unsaturated state caused by the release of confining pressure during the sampling process, samples were secondarily saturated using the vacuum saturation method for more than 72 h. The porosity of selected samples was determined using a cutting ring with a specific volume (119.9 × 10−6 m3). Eight samples from the borehole SY1 and 25 samples from the borehole SY2 were dried at 105°C for 48–72 h for analysis of their total porosities (), as shown in Figure 4.

3.2.2. Hydraulic Conductivity

The values of the samples in the profiles were determined using a laboratory testing method on undisturbed samples. A laboratory hydraulic conductivity tester was designed on the basis of a TST-55 permeameter, consisting of a cutting ring, two porous plates, a lantern ring, a top cover, a bottom cover, several screws, and pressure-resistant water supply bottle [38]. An external vacuum compressor was added to ensure that the of low-permeability clay was measured with high efficiency (Figure 5). The tester functioned within the operating pressure range of 0–0.6 MPa which was controlled with a precision regulator with an analytical precision of 0.02 MPa. The water flow in the low-permeability samples approximately followed at a constant pressure head. The external pressure generated by the air compressor was the most influential factor controlling water flow in the laboratory experiment. The change in the water supply bottle level was much lower than the external pressure head, and its function could be neglected (verifying test has been done but not present in this paper). Hydraulic pressure was produced in the water outlet if water flowed through the samples quickly during the experimental process, and then drip leakage was allowed in the test. For high-permeability samples (e.g., sand or silt deposits), an external vacuum compressor was ineffective, and the pressure head was provided by the water level of the supply bottle. was calculated using falling down head equations in accordance with highway engineering test methods for soils (2007) [39]. Testing results show that of the clay-rich zone varied with depth in the borehole SY1 from 7.3 × 10−10 to 1.5 × 10−11 m/s, and of silt-rich samples were between 2.6 × 10−8 and 1.8 × 10−9 m/s (Figure 4).

3.2.3. Diffusion Testing

Many different laboratory diffusion testing techniques for low-permeability materials have been reviewed [40]. The radial diffusion method was selected to apply in this study, because of its suitability for low-permeability aquitards and its high efficiency [4143]. The radial diffusion device was developed with special stainless steel, consisting of a tube, the upper cap, and the lower cap [44]. The height and internal diameter of the tube were 10 cm and 7 cm, respectively, as shown in Figure 6. A hole (length of 9 cm and diameter of 3 cm) was drilled along the central axis of the sample as a reservoir. Subsequently, water with a certain Cl concentration was injected into the reservoir. The application of Cl online monitoring in the reservoir assisted the recording of the concentration variation. To ensure that the probe of Cl online monitoring contacted the water in the reservoir easily, a spiral hole with the diameter of 3.0 cm was drilled on the upper cap.

Diffusing testing results were analyzed using the COMSOL software (COMSOL 4.4, Burlington, MA, USA) [45, 46] according to the initial and boundary conditions of the experimental model. The measured of Cl was 4.5 × 10−10 m2/s for a SY1 core sample (porosity, 0.4, depth 54.95–55.15 m, and laboratory temperature 12°C), and 3.5 × 10−10 m2/s for a SY2 core sample (porosity, 0.5, depth 32.8–33.0 m, and laboratory temperature 24°C). The measured values were corrected to the mean groundwater temperature (14°C) on the basis of the relationship between temperature and viscosity [47]. The corrected values of Cl for borehole SY1 and SY2 samples were 4.3 × 10−10 m2/s and 2.7 × 10−10, respectively.

4. Results and Discussion

4.1. Salt Source of Porewater

Salt in groundwater derives from several possible sources. The relationship of Cl and Br was always applied to determine the origin of groundwater salinity, since the Cl/Br ratios keep constant during the dilution and evaporation process of seawater prior to halite precipitation [48, 49]. Besides, Cl concentrations together with δ18O and δ2H values were further to characterize the salinization process of groundwater [50].

The Cl/Br mass ratios of SY1 saline porewater have a mean value of 203.6 below the near-shore seawater values. For SY2 saline porewater, Cl/Br mass ratios scatter around the near-shore seawater values with a mean value of 315.7. The quasilinear relationship of Cl and Br (Figure 7(a)) and the plot of Cl and δ2H (Figure 7(b)) suggest saline porewater salinity was derived from seawater. Most of the saline porewater plotted closely to the diluted line of standard seawater, indicating that the original porewater has been mixed by freshwater. The relationship of Cl and δ2H in SY1 suggests that the low Cl/Br mass ratios were probably attributed to the release of Br from diagenesis of marine organic material [51]. Thus, Cl could be assumed to behave conservatively and serve as a good natural tracer to quantify mechanism of solute transport.

4.2. Chloride–Depth Profiles

The extracted porewater Cl distributions of two boreholes along with the profiles were presented in Figure 4. The Cl concentrations yielded a well-defined 1D depth profile and indicated the presence of three distinct hydrogeologic zones interrupted by the aquifers. Two gradual decreases in Cl concentrations of borehole SY1 were observed from 334.8 to 60.0 mg/L (depth 162.2–86.5 m in zone B) and 16,086.1 to 262.2 mg/L (depth 19.8–2.0 m in zone A) as the depth decreases. Between 62.2 and 19.8 m, Cl concentrations increased rapidly from 16,086.1 to 412.3 mg/L. In borehole SY2, the high Cl concentrations were presented steadily ranging from 13,000 to 15,576.6 mg/L (mean 14,223.8 mg/L, depth 55.1–10.6 m). Then Cl concentration decreased with depth and attained its minimum value of around 250 mg/L at 85.6 to 120 m. The bad agreement between Cl concentrations in SY1 and borehole SY2 suggests that the chemical vertical variability in the profiles was assumed to be complex and large over distances [19].

The decreasing trend in observed Cl concentrations in the zone A from 19.8 m to the top of borehole SY1 could be attributed to the diluted mixture of modern meteoric water after the Holocene transgression. The shape of the Cl profile through the remaining zone A and the underlying deep first confined aquifer (aquifer I, Figure 4) suggests that Cl diffused from the zone A downward into the aquifer and upward toward the top of zone A. The lower Cl values (about 60 mg/L) from 108 to 89 m, attached to aquifer I, imply that the invasion of the different glacial meltwater into the aquifer I resulted in porewater in the underlying aquitards (zones B and C) suffering from vertical solute diffusive mixing.

The Cl concentrations of zone A in borehole SY2 were centered on the sand–rich deposits, especially for the depth of 54.2–10.6 m. It yielded a maximum Cl concentration of 15,576.6 mg/L, similar to the Cl concentrations of near-shore seawater values (about 9,000–16,000 mg/L). This could explain that the Holocene seawater rapidly travelled though the deposits and replaced the original porewater due to the coarse logic section. The B zone yielded Cl concentrations that were relatively constant at around 250 mg/L and showed a signature of freshwater, consisting with the constant Cl concentrations in borehole SY1 C zone and aquifer II (Figure 4). It shows that the deep first confined aquifer did not exert a significant impact on the adjoining Cl aquitard profile in borehole SY2 due to its thin thickness (about 3 m, depth 85.6–82.6 m). Due to their freshwater signatures, porewaters in the B and C zone of the two boreholes were assumed to derive from postdepositional glacial meltwater recharge with different time in the Late Pleistocene relative to groundwater in aquifer I.

4.3. Conceptual Transport Model

Vertical solute transport in a saturated medium can be simulated using the advection-diffusion equation: where is the average linear porewater velocity, D is the coefficient of hydrodynamic dispersion, C is the mass concentration of the solute, Z is the distance, and is the time. The coefficient of hydrodynamic dispersion is defined as D = + , where is the effective diffusion coefficient and α is dispersivity. The average linear porewater velocity is defined as V = , where is the Darcy velocity and is the effective porosity. We assumed that advection in the aquitard-aquifer system follows Darcy’s law and no threshold hydraulic gradient exists in the porewater transport process [52]. The effective porosity was supposed to be equal to the porosity, namely, (total porosity) [5, 42]. More, compaction-driven flow was also ignored according to the previous researches on aquitard-aquifer system [14, 15, 53]. Beyond that, the aquifer and aquitard were assumed to be individually homogeneous.

The porosities (s) for the borehole SY1 are uniform with depth (except depth of 14.82–15.02 m) and have a mean value of 0.41. For borehole SY2, s of aquitards are mostly around 0.5 and approximately 0.4 for the aquifers. Then, s of the aquitards and aquifers of borehole SY1 are set to be 0.41. Correspondingly, 0.5 and 0.4 are set to be s of the aquitards and aquifers in borehole SY2, respectively.

Based on the estimated values of two borehole aquitards using laboratory diffusion testing technique, the corrected values of Cl for boreholes SY1 and SY2 samples are presented in Table 1. The values for Cl in the aquitards and aquifers within 162.2 m were assumed to be 4.3 × 10–10 m2/s since no indication of pore structure occurred in these units [6]. The estimated for borehole SY2 was assumed to be 2.7 × 10–10 m2/s.

4.4. Initial and Boundary Conditions

The vertical distributions of Cl in the aquitard-aquifer system at the time of their deposition are unknown. The Cl profiles are probably attributed to several historical palaeohydrological events. The higher Cl concentrations in zone A originated from seawater during the Holocene transgression, and the low concentrations in aquifer I of SY1 possibly resulted from the invasion of glacial meltwater. However, the time of this intrusion of melt water is uncertain. The intrusion also had been described as aquifer “activation” of Hendry et al. (2013) [6] and Hendry and Harrington (2014) [11]. Defining the time of aquifer I recharge needed to depend on a great of assumptions. It is supposed that the recharge time for aquifer I was selected to coincide with the onset of the last glacial maximum, and the simulations were developed from 25 to 15 ka BP [29]. The relative constant values (around 250 mg/L) at depth of 159–120 m in borehole SY1 indicate that the aquitards and aquifers at depth of 182–19.8 m were assumed to be filled with freshwater (Cl, 250 mg/L) before “activation” of aquifer I. Afterwards, the later invasion with glacial meltwater (the last glacial maximum) into aquifer I resulted in Cl diffusion mixing transport at the top of zone B and the bottom of zone A. In fact, no realistic physical mechanism, other than perpetual flushing of the aquifer could be presented by a constant value throughout the aquifer I. A near-instantaneous change in Cl concentration within the aquifer I was assumed during the activation of aquifer [6]. And an instantaneous value of 60 mg/L was used in aquifer I after the activation of aquifer in keeping with the lower porewater values in the adjacent aquitard and groundwater values in aquifer I of the study area [20].

The sea level of the last glacial maximum was −130 m below the modern coastline [54]. In the Holocene, the sea level rose rapidly, and the study area was submerged by seawater [55], which seemed to result in the porewater salinization in the deposits, as Yi et al. (2012) [31] described occurring in the North China Plain. Due to the uncertainty of sedimentation rate and the relative thin thickness of the Holocene deposits, the fixed upper BC was set to locate at the surface of the Late Pleistocene since the onset of the Holocene, which was assumed to be the sediment-water interface (19.9 m). Theoretically, the Cl concentration for the upper BC should be the value of standard seawater (19.0 g/L). Actually, most of the near-shore seawater had low Cl concentrations (approximately 10,000–16,000 mg/L [20]), and the highest Cl concentration of 16,086.1 mg/L was observed in the two boreholes. Porewater salinity distribution in coastal plain aquitard-aquifer system was probably complex and mostly different from the standard seawater [15]. The assumed fixed Cl concentration of 16,100 mg/L was adopted as the upper BC.

The study area stayed in a marine environment until the appearance of Yellow River captured Huai River in AD 1,128. Since then, seawater retreated from the study area, and the shoreline gradually reached the present level (Figure 1(b)) [26]. It could believe that the Holocene deposits had been filled with seawater during the period of transgression, and the Cl concentration of seawater was set to be 16,100 mg/L. A freshening upper BC (fixed Cl concentration of 100 mg/L) was applied to the oxidation-nonoxidation interface (about 2 m) since AD 1,128, coinciding with the time of marine regression of this area.

For the borehole SY2, the depth of this study extends to the underlying aquifer III at the depth of approximately 160–180 m on the basis of the distributions of regional aquifers (Figures 2 and 4), because the drilled borehole did not reveal the underlying confined aquifer which should serve as the lower BC. Before the “activation” of aquifer I, the aquitards and aquifers at depth of about 160–18.9 m were assumed to be filled with freshwater (Cl, 250 mg/L). Subsequently, the glacial meltwater (the last glacial maximum; Cl, 60 mg/L) and the Holocene seawater intrusion led to Cl diffusion mixing transport of this profile. At the beginning of the Holocene, the initial Cl concentration of porewater at depth of 54.2–18.9 m was assigned as 16,100 mg/L since the relative coarser sediments occurred in this profile (described in Section 4.2). Beyond that, the similar various initial conditions and BC with different timeframes were assigned consisting with the setting for borehole SY1. Specific initial conditions and BC could refer to Table 1 for a summary of the boundary conditions.

4.5. Simulations of Aquitard-Aquifer System Cl Profiles

Previous researches on solute transport in the thick, regional aquitards have suggested that thin, sand-filled layers or permeable conduits might result in facilitating Cl “halos” laterally through an aquitard [56] and major deviations from one-dimensional diffusion profiles [12]. In this study, the sand layers (aquifers or sand streaks) are commonly encountered in the aquitards. Determining the influence of sand layers on the Cl transport simulations is crucial to palaeohydrological interpretation of porewaters.

The solute transport mechanism could be inferred to be advection if solute concentrations had no obvious vertical change; otherwise, curvilinear concentration profiles were attributed to molar diffusion transport [57, 58]. In order to detect the influence of sand layers on Cl transport, vertical 1D diffusion transport (see (1)) was postulated to be the dominant migration mechanism throughout the aquitard-aquifer system, because obvious vertical distributions appeared in the Cl concentration profiles (Figure 4).

As shown in Figure 4, the thickness of aquifer has an obvious influence on the Cl concentration profiles, for example, aquifer I (86.5–62.6 m) of SY1 borehole. The thick aquifer I is observed in borehole SY1, and Cl concentrations in the neighboring aquitards display the notable diffusion mixing trends occurring in the aquitard-aquifer interfaces. However, when the aquifer is relatively thin (the deep first confined aquifer with the thickness of about 3 m in the borehole SY2, depth of 85.6–82.6 m), the impact on Cl concentration profile in the neighboring aquitards is so little that even can be ignored. Namely, the thick aquifer should be regarded as the fixed concentration BC since the diffusion front does not reach the bottom of the aquifer in the simulations, while, for thin aquifer, the diffusion front reaches the bottom of aquifer in a short time, and the corresponding Cl concentration in the thin entire aquifer will keep changing with time on the basic of the Cl concentrations of the adjacent aquitards [15]. Besides that, the relatively thick aquifer (silt) is also found at depth of 52.3–41.5 m in borehole SY1. Groundwater in the aquifer entrapped palaeoseawater of the Holocene was similar to the constituents of the adjacent aquitards and was renewed more rapidly than the deep confined aquifers due to their shallow depth [19]. The Cl concentrations in the aquifer are also dynamic just like values in the thin sand layers. Thus, in this study, the thin aquifer II, the shallow first confined aquifer in borehole SY1, and the deep first confined aquifer (mentioned above, 85.6–82.6 m) in the borehole SY2 are supposed to be regarded as aquitards in the simulations.

To determine whether palaeohydrological change events affect the interpretation of the measured Cl profiles, the various boundary conditions were evaluated and applied to the modeling simulations. The transient constant upper and lower BC, initial conditions, aquifer boundary encountered in aquitards, and ’s were shown in Table 1. The starting time of aquifer I was reconstructed broadly varying from 30 to 12 ka BP to assess the impact of the last glacial maximum on the measured profile.

As is illustrated in Figure 8(a), the simulations with 20–15 ka BP are the best fit to the measured data at the depth of 162.2–62.6 m in borehole SY1. Using 15 ka BP simulation as the new initial condition and applying a subsequent saline BC (16,100 mg/L) with the transport time of about 10 ka and 70 ka are to develop the historical SY1 Cl concentration profile in the Holocene (Figure 8(b)). The obtained best fit simulation with transport time of 10 ka is consistent with the onset timescale of the Holocene. Subsequently, utilizing fresher upper BC with Cl concentration of 100 mg/L (the mean concentration of rainfall) which is referenced from Zhang et al. (2000) [59] and Zhang et al. (2003) [60] followed by saline BC, this model is used to simulate the shorter timescale (since AD 1,128). It provides a very good fit to the measured porewater Cl profile (Figure 8(c)), suggesting that the observed relatively large vertical scatter in Cl profile is closely related to the various upper BC.

Due to the thin deep first confined aquifer existing in borehole SY2, this aquifer activation is not considered in the simulations (note: the simulations related to the deep aquifer activation are modeled but not presented). As Figure 8(d) shows, the meltwater intrusion into the aquifer I has less of an impact on the Cl profile. According to the simulative result of the SY1 (Figure 8(a)), 15 ka BP simulation was selected as the new initial condition, and a series of simulations were performed, first by implementing a saline phase in the surface of the Late Pleistocene and then a freshening phase (Cl, 100 mg/L) in the oxidation-nonoxidation interface (about 2 m). Apart from the low Cl concentration in aquifer I (depth of 85.6–82.6 m), good fits between the simulated and measured Cl values are yielded (Figures 8(e) and 8(f)). These results could be explained that the realistic assumption is traceable in reasonable selection of the BC and initial conditions associations with the palaeohydrology events in the simulations.

4.6. Effects on the Diffusion Simulations of Cl Profiles
4.6.1. Effects of Effective Diffusion Coefficient on the Evolution of the Cl Profiles

Although simulations have the capacity to reflect tracer transport time, the simulated profiles cannot be unique because of the selection of values to material properties [61]. Sensitivity analyses of parameters in solute transport processes in aquitards have been conducted, and the results indicate that the diffusion model was most sensitive to the choice of values [62]. However, obtaining values using field or laboratory test methods is difficult because of the heterogeneity of deposits, test operations, and sampling. The selected unsuitable parameters for the simulation possibly lead to deviations in transport time [10]. In this study, the values of the sand layers (aquifers and sand streaks) are estimated similar to the aquitards in the two boreholes due to the absence of laboratory measurements of sand layers. To illustrate the influence of of the aquitard-aquifer systems on the simulated Cl concentration profiles, different values are selected for aquitard and aquifer deposits. In order to describe convenience, the simulated results of SY1 (Figure 8(c)) and SY2 (Figure 8(f)) are called “base case simulations,” as illustrated in Figure 9.

As the results shown in Figures 9(a) and 9(c), there is no greater difference between the base case simulation and the reconstruction simulations due to the thin sand layers in the borehole SY2. These also exhibit that same for the aquitards and aquifers would not result in the obvious offsets of the Cl simulation profiles in this study. The obvious differences were observed in the Figures 9(b) and 9(d), suggesting that the Cl simulations in the aquitard-aquifer system are affected by different chosen values. The best fits are obtained when = 4.3 × 10–10 m2/s for SY1 and = 2.7 × 10–10 m2/s for SY2 core samples, coinciding with the measured values.

4.6.2. Effects of Velocity Coefficient on the Evolution of the Cl Profiles

Advection transport is another factor that may influence solution migration [3, 62]. On the basis of the dynamic data of groundwater table monitored in the study area for numerous years, the mean water table level in the study area is generally 0.4 m; the mean piezometric level of the first and second confined aquifers is −4.5 m and −15.4 m, respectively. The determined vertical hydraulic gradient is fairly uniform with a value of 0.25 across the aquitards, assuming Darcy’s law is valid at very low velocities. Accordingly, the average downward porewater velocities through the clay-rich deposits are calculated to be between 0.118 m/ka and 5.76 m/ka for the values determined through the lab experiment method. In these simulations, groundwater velocities are supposed to range broadly from 0.01 to 10 m/ka.

The simulated profiles for downward velocities of 0.01, 0.1, 0.5, 1, and 10 m/ka and the measured Cl concentrations of the two boreholes are presented in Figure 10. The best fits simulated profiles for the two boreholes are obtained for velocity of 0.1 m/ka and less. The higher velocity cases evidently depart from the measured Cl profiles. However, the added advective velocity produces no marked change in the simulated profiles, and the minimum deviation does not affect the transport profiles. Therefore, it can surmise that diffusion solely is adequate to reconstruct the transport of Cl in the aquitard-aquifer system.

However, the Cl concentrations of porewater in the aquitards which are below aquifer I, and above aquifer III in boreholes present significant horizontal distribution and they do not show a coherent diffusion trend, it is likely that the strong overpressure in aquifers results in an advective component to the reversal of the Cl profiles [3]. Besides that, lateral advection was also not involved in these simulations, which probably contributes a horizontal advective component to the Cl profiles [63].

5. Conclusion

Numerical simulations of diffusive Cl transport were used to define the recharge processes of the aquifer I and explore the long–term evolution of porewater chemistry in aquitard-aquifer system of NJCP. Based on the preexisting knowledge of the sedimentary depositional environments and the events of transgression and regression since the Late Pleistocene, the various initial and upper boundary conditions over geological time were implemented. Diffusion has been approved to be the dominant transport mechanism of Cl in the aquitard-aquifer system. The thickness of the sand layers had significant influence on the shapes of the simulated Cl concentration in the aquitards. The simulations yielded well-defined 1D Cl profiles, suggesting that the assumed geological timeframes for the simulations are applicable and the high Cl profiles through the upper aquitards of the two boreholes are consistent with the onset of the Holocene (10 ka BP).

The results of this study present an understanding of the vertical solute transport in the aquitard-aquifer system based on the measured tracer profiles and transport parameters. The good agreement obtained between simulated and measured profiles suggests that the simulations can not only help one to understand solute transport mechanism in the aquitard-aquifer system but also provide insight into the timing of major geologic events (e.g., glaciations, transgression). The model simulation mechanisms are also suitable to other areas of the new Silk Road with the similar hydrogeological characteristics of NJCP.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This study was funded by the National Natural Science Foundation of China (Grants nos. 41502231 and 41272258). The authors would like to thank Professor Zhang Wen and two anonymous reviewers for their helpful comments and advice, as well as all the members of the Geological Survey of Jiangsu Province, China, for their assistance and support in the process of sampling and data collection.