#### Abstract

The dynamic behavior of groundwater flow and salt transport is affected by tide and pumping in coastal multilayered aquifers. In this paper, two groups of experiments were conducted considering different constant head inland boundaries. The fluctuation of the groundwater level and the process of seawater intrusion in the multilayered aquifers were observed. A two-dimensional SEAWAT model is developed to simulate the seawater intrusion to coastal aquifers under the influences of tidal fluctuation and groundwater exploitation. The hydrogeological parameters in the model are calibrated by the records of the groundwater level and salinity measurements. The results showed that the simulated groundwater level and salt concentration match the observation well. The groundwater level has the characteristics of periodic fluctuation with tide. The lag time of the groundwater level fluctuation in each monitoring point increases slightly with the increasing distance from the saltwater chamber. For the low tide, the inland freshwater recharge has main effect on groundwater level fluctuation. The rising tide has a negative effect on the drawdown of the groundwater level induced by pumping. For the high tide, the tide plays a major role on groundwater level fluctuation, compared with the inland freshwater recharge. Compared with the condition of high head of inland recharge, larger saltwater intrusion lengths and area have been observed and simulated in the aquifer, which means that faster inland motion of the seawater wedge would occur when the inland recharge is small in the coastal aquifers. It revealed that inland recharge plays a major role in the seawater intrusion for the same pumping rate of groundwater in different seasons. The analysis provides insights into how the tide fluctuation, groundwater pumping, and inland recharge effect on the area and rates of seawater intrusion.

#### 1. Introduction

Coastal areas are the most active area of human economic and social activities, which is built by many cities and large projects [1, 2]. The physical and chemical equilibrium of the coastal aquifer system is easily destroyed by the human activities, resulting in the problems of land subsidence, seawater intrusion, and environment deterioration. Among them, the seawater intrusion is a problem of worldwide concern, which is caused by the overpumping of groundwater, sea level rising, change of climate, and change of land use in the coastal area [3–9]. Among these factors, groundwater pumping is considered to be one of the most important challenges that promote the extent and severity of seawater intrusion [6, 9, 10]. As a result, the groundwater salinization can be directly caused by seawater intrusion due to the exploitation of groundwater [11].

The extent of seawater intrusion has been largely aggravated by the overexploitation of groundwater in coastal areas around the world, especially in China [6, 12–15]. For example, overexploitation has resulted in seawater intrusion in the Laizhou Bay, China [16]. The salinities of many groundwater samples were greater than 1000 mg/L up to tens of kilometers inland from the coastline, and the area affected by seawater intrusion has reached larger than 700 square kilometers [17, 18]. Additionally, in the Jiaozhou Bay area of China, the seawater intrusion has brought great damage in local agriculture and industry production [14]. Therefore, the mechanism, prediction, and prevention of seawater intrusion should be studied in the coastal zone.

As is known to us, the seawater intrusion process is a complex problem, because the seawater intrusion mixing with fresh groundwater couples water flow and salt transport. It involves the variable density groundwater flow, which is difficult to simulate well. Two types of flow models were developed for seawater intrusion simulation, including the sharp interface models and mixing interface models [14, 19–24]. There are many analytical solutions for variable density flow models to solve the coupled seawater intrusion equations under special conditions [6, 19, 20, 25]. Most of the analytical solutions are predominantly based on the Ghyben-Herzberg relation in the sharp interface models, in which the freshwater and saltwater are treated as separated fluids by an interface. As a consequence, the models may overestimate the actual penetration length of the seawater wedge [26]. Therefore, the analytical models of variable density flow including mixing between freshwater and seawater were provided considering the realistic condition [24, 27]. However, the analytical solutions are often derived based on some hypothetical conditions. The numerical approaches have been widely used to simulate seawater intrusion in the mixing interface [28, 29], which includes FEFLOW [30], SUTRA [31], and SEAWAT [14, 32, 33]. However, most of the numerical simulations were developed in two-dimensional cross sections, owing to the limitation of computation.

The dynamic behavior of groundwater flow induced by tidal fluctuation is an important issue in nature. There are some analytical solutions and numerical models adopted to study the influence of tides on the groundwater fluctuations and seawater interface [34–39]. The results show that the groundwater level fluctuates with a periodicity similar to tide periodicity, and the fresh saline water interface itself fluctuates due to sea level oscillations. However, the dynamic behavior of groundwater flow in a coastal aquifer is complex when considering the effect of groundwater exploitation [12, 16, 40–42]. Due to overpumping groundwater in coastal areas, the depression cones of the groundwater level can exacerbate seawater intrusion by reducing the hydraulic pressure of the freshwater. In addition, the inland boundary conditions have effects on seawater intrusion in coastal aquifers [43–46]. However, the comprehensive influence of tidal fluctuations, groundwater exploitation, and inland recharge on the location of the fresh saline water interface, including the groundwater flow in coastal aquifers, has received less attention. Moreover, the coastal aquifers are a multilayered aquifer system in most cases, and there are few studies on the seawater intrusion in the multilayered aquifers.

The objective of this study is to investigate the effects of a fluctuating sea level, groundwater pumping, and inland recharge on seawater intrusion in coastal multilayered aquifers. In this paper, two groups of experiments were conducted considering different constant head inland boundaries. The fluctuation of the groundwater level and the process of seawater intrusion in the unconfined aquifer and confined aquifer were observed. The numerical results were analyzed and compared with the experiment data. The effects of aquifer heterogeneity, boundary condition, groundwater pumping, and tide on solute mixing of saltwater and freshwater are described. The range of the saltwater and freshwater mixing zone and seawater wedge toe motion in the unconfined aquifer and confined aquifer are determined. The hydrogeological factors affecting salt distribution in the multilayered aquifer system are identified.

#### 2. Experimental Setup and Test Cases

##### 2.1. Experimental Setup

Figure 1 shows the laboratory tank, which is 6.6 m long, 0.6 m wide, and 1.5 m high, made of plexiglass with the thickness of 1.5 cm. The tank consists of three parts, which are the saltwater chamber, medium area, and freshwater chamber along the length direction. The medium area and the side chambers at either side are separated by two filter plates with small holes. The filter plates on both sides are 30 cm away from the left saltwater side and the right freshwater side, respectively. Geotextiles are laid on both sides of the filter plate, which acts as a barrier to sand passage. In order to avoid the deformation of the tank after filling sand, the main body of the tank is protected by a steel frame in the interval of 1 m. The seawater area and freshwater area are connected with the reservoirs, respectively. The two reservoirs are made of PVC boards, with the length of 140 cm, width of 120 cm, and height of 60 cm, placed on the ground. The right freshwater chamber is connected to a fixed head device (an overflow device for controlling the water level), which makes the freshwater head a fixed value. The left seawater chamber is controlled by the power cabinet system (controlling the tide device). The power cabinet system can control the inlet pump or outlet pump, which is connected to the seawater area, to make the seawater level fluctuate tidally.

The experiments required the measurements of the water level and the concentration of salt. The sea level on the left side of the tank was measured using a meter stick fixed on the steel frame. There are twelve monitoring point of pressure measuring tubes placed at the bottom of the back of the sand tank to measure the water level variation with time. In order to prevent the fine sand particles from entering into the pressure measuring tubes, a filter with a mesh diameter of 0.18 mm is placed at the connection between the piezometric tube and sand tank; the interval between every two pressure tubes is 50 cm in the horizontal direction. In order to study the seawater intrusion under the effect of groundwater exploitation, two groups of wells are set in the middle of the tank, which are 4 m and 5 m far from the saltwater boundary. The shallow wells No. 1 and No. 3 are installed in the fine sand layer, with the depth of 55 cm. The bottom of the wells was placed at the interface between the fine sand layer and clay layer. The deep wells No. 2 and No. 4 are installed in the coarse sand layer, with the depth of 156 cm. The bottom of the wells was placed at the bottom of the tank. The two groups of wells are cylindrical pipes, made of PVC material, with the diameter of 80 mm. There is a silastic tube with 2.79 mm internal diameter in each deep well, which connected to a peristaltic pump head (YZ35-13 type, 3 rollers, speed less than 600 rpm, and max flow rate 11000 mL/min) to extract water at the desired flow rate. The PVC pipes in both coarse sand layer and fine sand layer were drilled with small holes. In order to prevent the fine sand particles or coarse sand particles from entering into the wells, the bottom and the holes on the sidewall of the wells are wrapped by a filter with the mesh diameter of 0.18 mm.

In order to track the variation of the length of the saltwater wedge and the migration of the seawater-freshwater interface, there are five multiport sampling wells installed in the sand tank, which are 0.6 m, 1.8 m, 3 m, 4.2 m, and 5.4 m far away from the saltwater boundary, respectively. The multiport sampling wells are made of PVC board and have 6 sampling ports. They are labeled ports from A to F from the bottom to the top; the interval distance between every two ports is 20 cm. Each port is connected to a tube, whose top was connected to a Luer lock three-way valve at the top of the well through a Typon tubing.

Considering the heterogeneity of a coastal aquifer, the main body of the tank is filled with fine sand, clay, and coarse sand from top to bottom, with the thickness of 40 cm, 25 cm, and 50 cm, respectively. The bank slope of the fine sand layer is set to 6°. The setting of the slope satisfied the requirement that it should be stable and should not collapse during the period of the experiment. The sand samples are taken from the field site, which is located in the coastal area of Longkou City (the longitude is 120°3142, and the latitude is 37°4423), China. The particle sizes of fine sand and coarse sand are ranging from 0.25 mm to 0.1 mm and from 0.5 mm to 1 mm, respectively. The particle sizes of the clay are less than 0.01 mm. The hydraulic conductivities of the fine sand and coarse sand are 4.32 m/d and 37 m/d, obtained from the slug test in the wells No. 1 and No. 2, respectively. The hydraulic conductivity of clay is 0.25 m/d, which is obtained from a variable head permeameter in the laboratory.

The bottom of the sand box is the datum level. The simulated average sea level is 0.85 m. The highest and the lowest values of the tide are 1.05 m and 0.65 m, respectively. The tidal range is 0.4 m, and the tidal cycle is 40 min. The deionized water is used to be freshwater, which is made in the laboratory. The concentration and density of the freshwater are 0 g/L and , respectively. The saltwater is configured by deionized water and sodium chloride (analytical purity). The concentration and density of the saltwater are 18 g/L and , respectively. The bright blue is dissolved with salt to trace the saltwater transport in the tank and monitor the seawater wedge visually, with the concentration of 2 g/L. During the experiment, the salinities were measured at various depths using the multiport wells. The shape of the seawater wedge is recorded every 10 minutes by taking photos.

##### 2.2. Test Schemes and Flow Chart

Two groups of experiments (Cases 1 and 2) were designed and conducted to investigate the effects of tide and groundwater exploitation on groundwater flow and salt transport in a coastal multilayered aquifer system. In order to investigate the impact of inland freshwater recharge on seawater intrusion, these two cases of constant inland boundary are considered. The constant head of inland boundary is 95 cm for Case 1 and 80 cm for Case 2, which represents the high and low groundwater levels in different seasons. There are three tests in each case. The experiment starts after the sand sample installation.

For each case, the tank was initially flushed using freshwater with the saltwater chamber closed. For Case 1, the water level of freshwater on the right side was set to 95 cm. The freshwater flowed from the right side to the left side of the tank, and the excess amount of freshwater discharged from the left chamber. When the freshwater level on the left side rises to 95 cm, the steady flow was formed on both sides. Then, the tide device was applied, and the freshwater in the seawater chamber was replaced quickly with dyed seawater. The sea level fluctuated periodically, with the amplitude of 20 cm and period of 40 min. The pressure heads at the twelve measuring holes were observed. The electrical conductivity of water samples at different depths was measured by the conductivity meter. When the shape and position of the seawater wedge do not change within 10 minutes, it is assumed that the interface between saltwater and freshwater has reached a steady state condition. Subsequently, the peristaltic pump in well No. 4 started to extract water at the flow rate of 1.0 L/min. The pumping experiment stopped until a new steady state of seawater wedge was achieved. At last, the groundwater was continuously pumped with a larger flow rate of 2 L/min. The experiment stopped until the seawater wedge reached a new steady state. During the whole experiment, the water samples in the monitoring holes were collected at the interval between 5 min and 30 min. The migration of seawater wedge was observed every 5 minutes by taking photos during this period.

The water level of freshwater on the right side was set equal to 80 cm for Case 2. The whole experiment process was similar to that of Case 1, which was not described in detail. During the experiment, the images of the saltwater intrusion were captured with a high-speed camera with a resolution of pixels and an 8-bit grayscale pixel depth. After the experiment, a MATLAB code was then used to analyze all the experimental images to calculate the toe length of the saltwater wedge and provide maps of the solute concentration in the aquifers.

#### 3. Numerical Simulations of Variable Density Flow and Solute Transport

##### 3.1. Mathematical Model

The mathematical model for groundwater flow with variable density is applied to simulate the seawater intrusion processes. Based on the finite difference groundwater flow simulation model MODFLOW, the numerical simulation model of solute transport in groundwater flow is established using SEAWAT, considering the effect of density on groundwater flow [32]. The software has been widely applied in seawater intrusion and submarine groundwater discharge (e.g., [12, 47–49]).

The governing equation for the variable density groundwater flow can be expressed as
where is the flow direction; are the hydraulic conductivities in different directions (LT^{-1}); is the specific storage in terms of freshwater head (L^{-1}); is the equivalent freshwater head (L); is the effective porosity of porous medium (-); is the density of saline groundwater at a point in aquifer (ML^{-3}); is the density of freshwater (kg/m^{3}); is density of water entering from a source or leaving through a sink (ML^{-3}); is the volumetric flow rate of sources or sinks per unit volume of aquifer (T^{-1}); is the solute concentration (ML^{-3}); and is the time (T).

The process of seawater intrusion is a variable density groundwater flow, including the seepage, dispersion, and diffusion. The partial differential term of the concentration changing with time is added in the flow equation (1), and the change of concentration affects the groundwater flow of the aquifer. Therefore, the seawater intrusion requires the coupling solution of groundwater flow equation and solute transport equation. The equation of solute transport includes the groundwater convection term, hydrodynamic dispersion term, source-sink term, and reaction term, which is expressed as
where is the hydrodynamic dispersion coefficient tensor (L^{2}T^{-1}), is the solute concentration of water entering from sources or sinks (ML^{-3}), is the average linear velocity (LT^{-1}), and is the reaction term of chemical substance.

##### 3.2. Numerical Simulation Method

The aquifer system was assumed to be heterogeneous and isotropic. The length, width, and height of the simulation zone were 600 cm, 60 cm, and 115 cm, respectively. This model domain was discretized into the hexahedron elements with . There were 3150 columns and 15 layers with 47250 cells in total. The top boundary was a free surface, and the bottom boundary was treated as a no-flow boundary. On the right side boundary, a constant hydraulic head of 95 cm for Case 1 (80 cm for Case 2) was set. The concentration of freshwater was 0 g/L. The flow boundary at the left side boundary was defined as a time series of the hydraulic head based on a simple harmonic tidal function:
where (L) is the initial saltwater hydraulic head, (L) is the tidal amplitude, and (Rad T^{-1}) is the tidal angular frequency. Salinity was prescribed to be 18 g/L on the sea floor. The simulation period was 1200 minutes. A constant time step of 1 min was defined for the simulations in order to enable the analysis for the results from the numerical model. For every time step, the hydraulic heads on the sea floor changed according to equation 3.

Initially, the transient simulation was used to determine the extent of seawater intrusion under the condition of tidal fluctuation, without considering the pumping in the aquifer. Then, a pumping well was introduced, and the simulation continued until the new steady state conditions were obtained in different pumping rates.

##### 3.3. Model Calibration

In our model, the parameters of aquifers were calibrated by a trial-error method repeatedly to fit the head and salinity measurements for Case 1. The process of calibration was to adjust the values of the parameters, including the hydraulic conductivities, specific yield, specific storage, effective porosity, and dispersivity, until a good agreement between simulated and observed results is reached. Based on the empirical data, the horizontal hydraulic conductivity was assumed to be ten times of its vertical value of the aquifer system. The transverse dispersivity () was one-twentieth of the longitudinal dispersivity () [50, 51]. The longitudinal dispersivity and transverse dispersivity were determined mainly based on model calibration. The parameters used for the simulations were given in Table 1. From the estimated values listed in Table 1, one can see that the hydraulic conductivity of clay is smaller, compared with those of the fine sand and coarse sand. Therefore, the clay is regarded as the semipermeable layer. The fine sand layer and coarse sand layer are the unconfined aquifer and confined aquifer, respectively.

Figure 2 shows general similar groundwater dynamics between simulated and observed groundwater level data in the confined aquifer during the simulation period of 120 minutes, despite discrepancies around the peaks and troughs. The relative error (RE) and correlation coefficient () are used to quantitatively analyze the fitting results, which can reflect the fitting degree between the calculated values and measured values. The calculated results are shown in Table 2. The RE values between observed and calculated hydraulic heads at the monitoring points are ranging from 5.79% to 13.94%, and the correlation coefficient () is between 0.8265 and 0.9659, respectively. The results showed that most of the RE values are less than 10% and the values are larger than 0.9. It indicates that the estimation of hydrogeological parameters is reasonable and the reliability of the model is high.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 3 shows the observed values and the calculated values of salt concentrations at the sampling wells in the confined aquifer during the simulation period. From Figure 3, one can see that the observed and calculated values of Cl^{-} concentrations increased during the period of 0 to 480 minutes. Then, they became stable from 480 minutes to 1200 minutes. The Cl^{-} concentrations simulated by the model were consistent with the observed ones in the sampling wells No. 1 and No. 2. The concentration of saltwater in sampling well No. 2 is lower than 18 g/L, which may be due to the fact that the freshwater flows into the saltwater chamber during the experiment. These errors are small relative to the maximum concentration variations in the sampling wells. Generally, the simulated groundwater levels and salt concentrations match the observations well.

**(a)**

**(b)**

#### 4. Results

##### 4.1. Effect of the Tidal Level on the Groundwater Level

In order to study the fluctuation characteristics of the groundwater level caused by tidal effect in the aquifers, Case 1 is selected to analyze the variation of the tide and the groundwater level of the confined aquifer, without considering the pumping effect. The period of tidal fluctuation was 40 minutes. The amplitude of the tide was 20 cm, and the average sea level was 85 cm. The monitoring period was from 0 min to 200 min, and the groundwater level fluctuated with time during the period of 200 min.

Figure 4 reports the variation of the groundwater level at the selected seven monitoring points and the tidal level fluctuation. The graph adopts the data from 0 minute for five consecutive tidal periods since the beginning of the experiment. From the figure, one can see that the groundwater level and tidal fluctuation have the same regularity, which shows the characteristics of periodic fluctuation. The farther the monitoring point is from the left boundary (saltwater chamber), the smaller the fluctuation range of the groundwater level is. There is a lag time between the peak or valley of the tide and that of the groundwater level. In addition, there is a certain difference in the time of the peak value or the valley value at each monitoring point. The lag time between the peak value of the groundwater level at each point and that of the sea level from the shoreline (monitoring point No. 1) to inland (monitoring point No. 12) is ranging from 5 min to 10 min. It indicates that the lag time of the groundwater level fluctuation in each monitoring point increases slightly with the increasing distance from the saltwater chamber.

##### 4.2. Effects of Pumping on the Tide-Induced Groundwater Level

In order to study the pumping effects on tide-induced groundwater level fluctuation in the coastal multilayered aquifers, two groups of experiments were conducted. During the two groups of experiments, the change of the interface between seawater and freshwater and the variation of the groundwater level at each monitoring point were observed until the system reached the final steady state. In this section, we only focus on the variation of the groundwater level.

Figure 5 shows the observed values and calculated values of the groundwater level at the low tide, average sea level, and high tide for Case 1. The match between the observed and calculated groundwater levels is well. As the tide rises from the low tide (65 cm) to high tide (105 cm), the groundwater level increases under the condition of tidal effect without pumping (blue lines of Figures 5(a)–5(c)). From Figure 5(a), one can see that the groundwater level decreases when the pumping rate equals to 1 L/min for the low tide. The groundwater depression cone is formed in the well No. 2 due to the pumping. The drawdown of the groundwater level in well No. 2 reaches the maximum with a value of 8.5 cm. The groundwater level drawdown on both sides of the pumping well No. 2 decreases gradually to the left boundary and right boundary. It indicates that transient pumping can significantly enhance the amplitude of the groundwater level fluctuation. A further increment of the pumping rate to caused faster drawdown of the groundwater level in the aquifer. The drawdown of the groundwater level in well No. 2 can reach the maximum value of 15 cm, compared with that of the case without considering pumping. From Figure 5(c), one can see that the groundwater level drawdowns in well No. 2 are 5.0 cm and 9.7 cm for the high tide, respectively, when the pumping rates are 1 L/min and 2 L/min. When the tide rises from low tide to high tide, the drawdown of the groundwater level at each point decreases (Figures 5(a)–5(c)). It indicates that the rising tide has a negative effect on the drawdown of the groundwater level induced by pumping.

**(a)**

**(b)**

**(c)**

In order to investigate the inland freshwater recharge on the groundwater dynamic, the freshwater head on the right boundary was set to 80 cm (Case 2), which is lower than that of Case 1 (95 cm). Figure 6 shows the observed values and calculated values of the groundwater level at low tide, average sea level, and high tide for Case 2. The simulated results captured the fluctuating trend of the groundwater level with time at the monitoring points. From Figure 6, it can be seen that the groundwater level increases under the condition of tidal effect without pumping, as the tide rises from the low tide to high tide. Thus, the variations of the groundwater level for Case 2 are similar to those of Case 1. It can be seen from Figure 6(a) that the groundwater level drawdowns in well No. 2 can reach the maximum values of 12 cm and 18 cm, respectively, when the pumping rates are 1 L/min and 2 L/min. The groundwater level drawdown on both sides of the pumping well No. 2 decreases gradually to the seaside boundary and inland boundary under the condition of different pumping rates. Compared with Figure 6(a), from Figures 6(b) and 6(c), one can conclude that the groundwater level drawdown decreases at the monitoring points, as the tide rises under different pumping rates. When the sea level is 105 cm, the groundwater level drawdowns in well No. 2 are 5.6 cm for and 8 cm for , respectively, which are smaller than those at low tide. From Figures 5(a) and 6(a), one can see that the groundwater level drawdowns at the monitoring points for Case 1 are lower than those of Case 2, when the pumping rates are 1 L/min and 2 L/min. Especially in well No. 2, the groundwater level drawdown for Case 1 is lower than that of Case 2, with values of 3.5 cm for and 3.0 cm for . It indicates that the inland freshwater recharge in Case 2 is smaller than that of Case 1, which cannot replenish groundwater exploitation in time. However, when the sea level reaches its maximum, the difference between the groundwater level drawdowns of Case 1 and those of Case 2 is small. In well No. 2, the drawdown for Case 1 is slightly smaller than that of Case 2 with the value of 0.6 cm for and larger than that of Case 2 with the value of 1.7 cm for . It reflects that the tide plays a major role in the high tide, compared with the inland freshwater recharge. A large amount of seawater enters into the aquifer, because the groundwater level is lower than the sea level.

**(a)**

**(b)**

**(c)**

##### 4.3. Saltwater Intrusion in the Multilayered Aquifers

During the experiments, there is no seawater intrusion in the middle semipermeable layer (clay layer), because of the low permeability. Therefore, in this part, we only focus on the analysis of the seawater intrusion processes in the unconfined aquifer (fine sand layer) and confined aquifer (coarse sand layer). The salinity at the mixing interface is defined as larger than 1% of salinity of seawater.

Figure 7 reports the transient experimental saltwater wedges at different time intervals until the system reached a steady state condition for Case 1. From Figure 7, it can be seen that the interface between saltwater and freshwater is not clear at the initial stage of seawater intrusion (about 80 min). The curved shape of the seawater wedges in the fine sand layer and coarse sand layer is linear. As the extent of seawater intrusion increases, the saltwater wedge migrates inland, and the interface between seawater and freshwater becomes clear. The curved shape of the seawater wedge changes from linear to a concave parabolic curve. The steady state condition occurred within 480 min (12 tidal cycles) in both models for Case 1 under the influence of the tidal level. There is an obvious interface between the saltwater and freshwater. Then, in order to investigate the effect of pumping on the seawater intrusion, the groundwater in well No. 2 was pumped out with the pumping rate of (1 L/min). The original equilibrium between the saltwater and freshwater was broken, and the saltwater wedge continued to move toward the land boundary until the secondary equilibrium was reached. Subsequently, the pumping rate in well No. 2 increases to (2 L/min) from the time of 760 min, and the shape of the seawater wedge remains a concave parabolic curve, until a new steady state was reached. During the process of seawater intrusion, the angles of seawater wedge tip in the fine sand layer and the coarse sand layer decreases gradually with time. The results show that the mixing degree of seawater and freshwater in the aquifer is small, and the transition zone between seawater and freshwater is narrow.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

Figure 8 presents the transient experimental seawater intrusion process at different time steps until the system reached a steady state condition for Case 2. Similar to Case 1, the interface between seawater and freshwater is not clear initially, and the curved shape of the seawater wedge is a straight line in the fine sand layer. As the tide fluctuates, the seawater wedge moves inland until it is stable, and the interface between seawater and freshwater becomes clear gradually. The shape of the seawater wedge changes from a straight line to a concave parabolic curve in the fine sand layer. The interface between seawater and freshwater is clear, which shows a shape of a concave parabolic curve in the coarse sand layer. The seawater intrusion in the unconfined aquifer and confined aquifer can reach a steady state at 500 min, which is relatively slower than that of Case 1. Then, in order to investigate the effects of pumping on the seawater intrusion, the groundwater in well No. 2 was pumped out with the pumping rate of (1 L/min) from 500 min and (2 L/min) from 840 min, respectively. During the pumping process, seawater intrudes toward inland continuously, and the shape of the seawater wedge still keeps a concave parabolic curve, until a new steady state was reached. The angles of seawater and freshwater interface tip in the fine sand layer, and the coarse sand layer decreases gradually. The widening of the seawater wedge for Case 2 was larger than that of Case 1. The reason is that the inland freshwater recharge is relatively small in Case 2 compared with that of Case 1.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

##### 4.4. Dynamics of Salt Wedge in the Multilayered Aquifers

Figure 9 shows the variation of salt wedge toe length and salt wedge area with time in the unconfined aquifer and confined aquifer for Case 1. Generally, the numerical model can depict reasonably well the development of the seawater wedge variation observed in the physical experiment. Time-varying toe length of Case 1 in the fine sand layer and coarse sand layer is depicted in Figure 9(a). From Figure 9(a), one can see that the transient toe length increases with time initially and then tends to be a stable value under the condition of tidal fluctuation. When the seawater intrusion reached the steady state, the transient toe length of the upper fine sand layer is 195 cm and that of the coarse sand layer is 242 cm (Table 3). The salt wedge migration rate toward inland decreases with time dramatically initially, and then, it decreases slowly to zero (Figure 9(b)). When the groundwater in well No. 2 was pumped with the rate of , the seawater further intrudes toward inland until a new steady state is reached. The transient toe length increases with time in the early state of pumping and then tends to be a stable value (Figure 9(a)). The intruding rate of seawater wedge increases with time and reaches a maximum value at the time of 540 min, and then, it decreases with time (Figure 9(b)). When the seawater intrusion reached the steady state, the transient toe length of the upper fine sand layer is 250 cm, and that of the coarse sand layer is 264 cm (Table 3). When the groundwater in well No. 2 increased to , the seawater further intrudes toward inland until a new steady state is reached. The transient toe length increases with time slightly until it is stable (Figure 9(a)). The salt wedge migration rate increases with time and reaches a maximum value at the time of 820 min, and then, it decreases to zero (Figure 9(b)). It indicates that the steady state of the seawater intrusion can be judged according to the migration rate of toe length. The maximum values of toe length are 261 cm in the upper fine sand layer and 272 cm in the lower coarse sand layer, respectively (Table 3). The toe length of the lower layer is larger than that of the upper layer, because the hydraulic conductivity of the lower layer is higher than that of the upper one. Figure 9(c) represents the time-varying curve of the seawater wedge area. Compared with Figure 9(a), one can see that the time varying of area is similar to that of length in the three stages. It indicates that the shape of the seawater wedge changes gradually during the process of seawater intrusion. The maximum area of seawater intrusion can reach 1.0 m^{2} in the fine sand layer and 1.36 m^{2} in the coarse sand layer, respectively.

**(a)**

**(b)**

**(c)**

Figure 10 presents the variation of toe length, migration rates, and area of the salt wedge with time in the unconfined aquifer and confined aquifer for Case 2. It shows that the numerical model predicted very well the inland motion of the saltwater wedge observed in the physical model. Similar to Figures 9(a) and 10(a) shows that the transient toe length increases with time initially and then tends to be a stable value under the condition of tidal fluctuation. The transient toe lengths of the steady state are 242 cm in the upper fine sand layer and 252 cm in the lower coarse sand layer, respectively (Table 3). Then, the transient toe length increases with time when the groundwater in well No. 2 was pumped with the rate of . When the seawater intrusion reached a new steady state, the transient toe length of the upper fine sand layer is 290 cm, and that of the coarse sand layer is 336 cm (Table 3). The transient toe length increases with time slightly until it is stable, when the groundwater in well No. 2 was increased to . Finally, the maximum values of toe length are 307 cm in the upper fine sand layer and 356 cm in the lower coarse sand layer, respectively (Table 3). The results show that substantially larger saltwater intrusion lengths have been observed and simulated in the unconfined aquifer and confined aquifer for Case 2 compared with those observed and simulated in Case 1. Figure 10(b) shows that the salt wedge migration rate toward inland decreases with time dramatically initially, and then, it decreases slowly to zero under the first stage. The seawater intruding rate increases with time under the condition of pumping and reaches a maximum value at the time of 560 min, and then, it decreases with time. Subsequently, it increases to be a maximum value at the time of 900 min, and then, it decreases to zero. Figure 10(b) also shows that the intruding migration rate in Case 2 was relatively larger than that in Case 1 in the three stages. This means that faster seaward motion of the seawater wedge would occur when the inland recharge is small in the coastal aquifers. Compared with Figure 10(a), from Figure 10(c), one can see that the time-varying curve of the area is similar to that of length in the three stages. The maximum area of seawater intrusion can reach 1.2 m^{2} in the fine sand layer and 1.78 m^{2} in the coarse sand layer, respectively, which are relatively larger than those in Case 1 during the process of seawater intrusion.

**(a)**

**(b)**

**(c)**

#### 5. Discussions

##### 5.1. Sensitivity Analysis

Different scenarios were investigated to assess the sensitivity of the seawater intrusion to the parameters’ values. The above-discussed simulation is referred to as the “base case,” taking the model of Case 1 as an example. The hydraulic conductivity and dispersivity are important parameters for seawater transport. Sensitivity analysis with respect to the hydraulic conductivity and dispersivity was conducted for the variation of the salt wedge. Only the values of the model parameter for sensitivity analysis were changed, and all the other parameters were fixed in each simulation.

In order to investigate the effect of hydraulic conductivity on the transport of salinity of seawater, the value of it for the unconfined aquifer and confined aquifer was changed from the base case. When the hydraulic conductivity was increased from 4.3 m/d (base case) to 43 m/d for the fine sand, and from 38 m/d (base case) to 380 m/d for coarse sand, the match to the observed toe length of seawater by the new model is worse than that of the base case (Figure 11(a)). The toe lengths simulated by the new model were consistently higher than the observed ones, but those of the base case were close to the observations in the upper and lower layers. The equilibrium time of the saltwater wedge of the new model is earlier than that of the basic model under the condition of tidal fluctuation and pumping. The intruding rate of the seawater wedge in the new model is larger than that of the basic model. The results indicated that the hydraulic conductivities of the two layers are high, and a large quantity of seawater flowed into the aquifer. When the hydraulic conductivity was decreased from 4.3 m/d (base case) to 0.43 m/d for the fine sand, and from 38 m/d (base case) to 3.8 m/d for coarse sand, the toe lengths of seawater simulated by the new model were obviously worse than those of the base case. From Figure 11(b), one can see that the arrival time of the steady state when decreasing the hydraulic conductivities of the two layers lagged behind that of the base case. Thus, one is led to conclude that the hydraulic conductivities of the two layers are too low, which are not able to reproduce the variation of the seawater wedge.

**(a)**

**(b)**

The sensitivity of toe length of seawater to dispersivity was explored to understand the contribution of dispersivity to salinization of the aquifer under pumping condition. In the simulations, the values of longitudinal dispersivity and the ratio of transverse dispersivity to longitudinal dispersivity were changed from the base case (, ). When the longitudinal dispersivity was increased to (keeping the dispersivity ratio constant), the match between the observed and simulated toe lengths became worse than those of the base model (Figure 12(a)). The toe lengths simulated by the new model were higher than those of the base model. However, the difference was not too large, suggesting that intermediate dispersivity values between the base model and new model would still be acceptable overall. When the longitudinal dispersivity was decreased to , Figure 12(b) shows that the toe lengths simulated by the new model were lower than those of the base model. The results showed that changing the dispersivity has no effect on the arrival time of the steady state saltwater wedge.

**(a)**

**(b)**

##### 5.2. Physical Mechanism

This study highlights the importance of the fluctuating sea level and inland recharge for pumping-induced groundwater salinization. The sea tide and inland freshwater recharge have impacts on the groundwater level and seawater intrusion process in the multilayered aquifers. The groundwater level has the characteristics of periodic fluctuation with the tide. The farther the position is from the sea, the smaller the fluctuation of the groundwater level is, which is similar to the results of previous studies [36, 49]. However, the groundwater level fluctuation is complex when considering the effect of pumping (e.g., [39, 42, 52, 53]). We should correct the drawdown data during the pumping test, due to the groundwater level dynamics jointly induced by tidal forcing and groundwater pumping. In contrast to those previous studies, this study concentrates on the comprehensive effects of the fluctuating sea level and inland recharge on the pumping-induced groundwater level. The groundwater level variation in the unconfined aquifer and confined aquifer was simulated, considering different constant head inland boundaries, which represents the high and low groundwater levels in different seasons. We found that the inland freshwater recharge has main effect on the groundwater level fluctuating under the condition of low tide. The drawdown of the groundwater level decreases as the tide rises, which indicates that the rising tide has a negative effect on the drawdown of the groundwater level induced by pumping. However, the tide plays a major role compared with the inland freshwater recharge when the tide reaches maximum. Our results show that the seasonal variation of inland recharge should be considered to improve understanding of basic controls on groundwater level variation due to pumping.

In addition, we considered the variation of seawater intrusion with time in different layers. The heterogeneity of the aquifer has great effect on the seawater intrusion process (e.g., [9, 19, 41]). Our results also show that aquifer heterogeneity has an effect on salinity distributions in the fine sand layer and coarse sand layer, due to the different hydraulic conductivities in the two layers. However, there is less attention on the combined effects of the sea tide and inland recharge on pumping-induced groundwater salinization. Our results provide some insights into the toe length, migration rates, and area of seawater intrusion under the condition of tide fluctuation, groundwater pumping, and inland recharge. We found that the inland recharge in different seasons plays a major role in the seawater intrusion for the same pumping rate of groundwater. When the constant head of inland recharge is large, there is much more freshwater flowing into the aquifer, which makes the saltwater restored. When the constant head of inland recharge is small, the length and area of seawater intrusion in the aquifer are relatively large, because they are governed by the hydraulic gradient between the sea level and groundwater level. The analysis provides insights into how the tide fluctuation, groundwater pumping, and inland recharge effect on the area and rates of seawater intrusion. For future work, further efforts to do more changes of the boundary condition effects on seawater intrusion should be considered.

#### 6. Conclusions

This study investigates the effects of a fluctuating sea level and inland recharge on pumping-induced groundwater salinization in coastal multilayered aquifers. A two-dimensional numerical model of variable density flow and solute transport is established. Using the laboratory experiments and numerical simulations, the fluctuation of the groundwater level and the process of seawater intrusion in the unconfined aquifer and confined aquifer were considered. The effects of aquifer heterogeneity, boundary condition, groundwater pumping, and tide on solute mixing of saltwater and freshwater are described. The main findings of the study are the following: (1)The farther the monitoring point is from the left boundary, the smaller the fluctuation range of the groundwater level is. There is a lag time between the peak or valley of the tide and that of the groundwater level. The lag time of the groundwater level fluctuation in each monitoring point increases slightly with the increasing distance from the saltwater chamber(2)In each simulation, the matches between the observed and calculated groundwater levels are well. The transient pumping can significantly enhance the amplitude of the groundwater level fluctuation. A further increment of the pumping rate caused faster drawdown of the groundwater level in the aquifer. For the low tide, the inland freshwater recharge has main effect on groundwater level fluctuation. The drawdown of the groundwater level at each point decreases as the tide rises, which indicates that the rising tide has a negative effect on the drawdown of the groundwater level induced by pumping. For the high tide, the groundwater level at the monitoring point is relatively higher when the head of inland recharge is large. It indicates that the tide plays a major role in the high tide, compared with the inland freshwater recharge(3)The transient toe length increases with time initially and then tends to be a stable value under the condition of tidal fluctuation and pumping for different inland recharges. Larger saltwater intrusion lengths and area have been observed and simulated in the unconfined aquifer and confined aquifer for Case 2 compared with those observed and simulated in Case 1. The intruding migration rate in Case 2 was relatively larger than that in Case 1 in the three stages. This means that faster seaward motion of the seawater wedge would occur when the inland recharge is small in the coastal aquifers. There is much more freshwater flowing into the aquifer for Case 1, which makes the saltwater restored. When the constant head of inland recharge is small, the length and area of seawater intrusion in the aquifer are relatively large, which is governed by the hydraulic gradient between the sea level and groundwater level. It revealed that inland recharge plays a major role in the seawater intrusion for the same pumping rate of groundwater in different seasons(4)The sensitivity analysis revealed that the estimated parameters of hydraulic conductivity and dispersivity are well determined. Increasing or decreasing the hydraulic conductivities of the two layers is not able to reproduce the variation of the seawater wedge. Changing the dispersivity has no effect on the arrival time of the steady state saltwater wedge. The analysis provides insights into how the tide fluctuation, groundwater pumping, and inland recharge effect on the area and rates of seawater intrusion

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

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

#### Authors’ Contributions

Qiaona Guo and Zhifang Zhou conceived and designed the study; Jiangwei Huang analyzed the laboratory data; and Qiaona Guo wrote the paper with the assistance of Jinguo Wang.

#### Acknowledgments

This research was supported by the National Key R&D Program of China (No. 2016YFC0402803), the National Natural Science Foundation of China (No. 41772235), and the Fundamental Research Funds for the Central Universities (No. 2019B16814).