Transport of Contamination under the Inﬂuence of Sea Level Rise in Coastal Heterogeneous Aquifer

: This paper provided for the ﬁrst time an experimental study on the inﬂuence of sea level rise on transport of contamination in the heterogeneous unconﬁned aquifer of the coastal zone. The experiments were conducted using the tank, considering the di ﬀ erence between sea level and inland head 1 cm for Case 1 and 2 cm for Case 2. Observed data were validated using the numerical model, which matched well with the toe length of seawater wedge and the shape of the contaminant plume. The results showed that the observed and simulated values of Cl − concentration at the sampling points increased sharply at the initial time, and then they increased slowly and tended to be stable. The seawater wedge migrated inland with time under the e ﬀ ects of the hydraulic gradient toward the inland and the density di ﬀ erence between saltwater and freshwater. The steady state length of the 50% isoline of the seawater wedge was 167 cm in Case 2, which was larger than that of Case 1. The maximum area of plume in Case 2 was 0.13 m 2 , larger than that in Case 1, which indicated that the velocity of di ﬀ usion of the contaminant plume increased as the sea level increased. As the velocity of di ﬀ usion increased, the time for pollutant migration to the intersection between seawater and freshwater became shorter. The maximum area and vertical depth of pollutant plume were sensitive to the hydraulic conductivity, dispersivity, and contamination concentration. The inﬁltration depth and range of the contaminant plume in the heterogeneous aquifer were greater than those in the homogeneous aquifer of the actual beach.


Introduction
The coastal area is densely populated, due to the rapid development of the economy there. Therefore, a large amount of groundwater is pumped to provide people with water for living and industrial production. Nevertheless, the problems of water resources, ecology, and environment pollution in coastal areas are complex and serious, due to the influences of interaction between land and sea, global warming, and sea level rising [1][2][3][4][5]. Among them, environmental pollution is a very serious problem. For example, the spilled oil deposits within the beach, caused by the accident of offshore oil tankers, which causes the long-term pollution of groundwater in the aquifer [6][7][8]. In addition, the wastewater infiltrates into the aquifer due to seepage from facilities, with the rapid development of the aquaculture in recent years, which results in the deterioration of groundwater quality and water balance of coastal aquifers [9][10][11]. Therefore, understanding the contaminate transport process and remediation in coastal groundwater is crucial to the economic, social, and environmental sustainability of coastal areas.
The spatial and temporal distribution of salinity in coastal zones have been focused on, because the seawater encroaches into the coastal aquifer e.g., [2,[12][13][14][15][16]. In the recent decade, there has been Sustainability 2020, 12, 9838 3 of 16 and seawater intrusion on contaminant transport are discussed. The main hydrogeological factors affecting the contaminant transport in the coastal aquifer are determined.

Laboratory Experiment
The dimension of the tank for experiment was 260 cm × 30 cm × 100 cm, which was shown in Figure 1. The tank was made of plexiglass and protected by a steel frame. It was composed of three parts: a freshwater chamber, a central chamber, and a seawater chamber. The left and right side chambers were separated from the central chamber by plates, which were made of PVC materials. In order to avoid the sand flowing from the central chamber to the two side chambers, the geotextile was pasted on the surface of the PVC board drilled with small holes. The bottom of the seawater chamber and freshwater chamber were respectively connected with a water reservoir, which was made of organic plastic plates with the length of 84 cm, width of 84 cm, and height of 84 cm. The seawater and freshwater were respectively pumped into the seawater chamber and freshwater chamber by peristaltic pumps. The left and right side chambers represented the freshwater and seawater boundary conditions, respectively. Three different particle sizes were used during the experiment, namely fine sand, medium sand, and coarse sand from up to down. The material was homogeneous in each layer. The particle sizes of fine sand, medium sand, and coarse sand were ranging from 0.1 to 0.25 mm, from 0.25 to 0.5 mm, and from 0.5 to 1 mm, respectively. The material in each layer was compacted. In reality, the coastal aquifer inclines to the sea boundary. However, the slope of the aquifer was generally small. Therefore, the slope of the aquifer for the experiment was set horizontal. The thickness of the surface sandy aquifer was generally larger, considering the actual beach condition. Therefore, the fine sand layer, medium sand layer, and coarse sand layer were laid horizontally in the sand tank, with the thickness of 23 cm, 16 cm, and 16 cm, respectively. The thickness of the whole aquifer was 55 cm. First of all, a certain amount of freshwater was put into the tank before filling the sand. Then, the coarse sand with the thickness of 2 cm was filled. The sand was left for 30 min to make it fully saturated. Subsequently, the above steps were repeated in turn until the coarse sand layer reach 16 cm. The method of laying sand in the medium sand layer and fine sand layer was similar to that of coarse sand. During the process of laying sand, each layer was compacted as much as possible to remove the air. The hydraulic conductivities of the fine sand, medium sand, and coarse sand were 2.3 × 10 −5 m/s, 2.9 × 10 −4 m/s, and 6.9 × 10 −4 m/s, respectively, which were obtained from the variable head permeameters in the laboratory.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 17 examined, based on the laboratory experiments and numerical simulation. The effects of boundary condition, aquifer heterogeneity, and seawater intrusion on contaminant transport are discussed. The main hydrogeological factors affecting the contaminant transport in the coastal aquifer are determined.

Laboratory Experiment
The dimension of the tank for experiment was 260 cm × 30 cm × 100 cm, which was shown in Figure 1. The tank was made of plexiglass and protected by a steel frame. It was composed of three parts: a freshwater chamber, a central chamber, and a seawater chamber. The left and right side chambers were separated from the central chamber by plates, which were made of PVC materials. In order to avoid the sand flowing from the central chamber to the two side chambers, the geotextile was pasted on the surface of the PVC board drilled with small holes. The bottom of the seawater chamber and freshwater chamber were respectively connected with a water reservoir, which was made of organic plastic plates with the length of 84 cm, width of 84 cm, and height of 84 cm. The seawater and freshwater were respectively pumped into the seawater chamber and freshwater chamber by peristaltic pumps. The left and right side chambers represented the freshwater and seawater boundary conditions, respectively. Three different particle sizes were used during the experiment, namely fine sand, medium sand, and coarse sand from up to down. The material was homogeneous in each layer. The particle sizes of fine sand, medium sand, and coarse sand were ranging from 0.1 to 0.25 mm, from 0.25 to 0.5 mm, and from 0.5 to 1 mm, respectively. The material in each layer was compacted. In reality, the coastal aquifer inclines to the sea boundary. However, the slope of the aquifer was generally small. Therefore, the slope of the aquifer for the experiment was set horizontal. The thickness of the surface sandy aquifer was generally larger, considering the actual beach condition. Therefore, the fine sand layer, medium sand layer, and coarse sand layer were laid horizontally in the sand tank, with the thickness of 23 cm, 16 cm, and 16 cm, respectively. The thickness of the whole aquifer was 55 cm. First of all, a certain amount of freshwater was put into the tank before filling the sand. Then, the coarse sand with the thickness of 2 cm was filled. The sand was left for 30 min to make it fully saturated. Subsequently, the above steps were repeated in turn until the coarse sand layer reach 16 cm. The method of laying sand in the medium sand layer and fine sand layer was similar to that of coarse sand. During the process of laying sand, each layer was compacted as much as possible to remove the air. The hydraulic conductivities of the fine sand, medium sand, and coarse sand were 2.3 × 10 −5 m/s, 2.9 × 10 −4 m/s, and 6.9 × 10 −4 m/s, respectively, which were obtained from the variable head permeameters in the laboratory.   The pollution source chamber was installed on the top of the aquifer, which was 106 cm away from the left freshwater chamber. It was made by the transparent organic plastic plates, with the length of 30 cm, width of 5 cm, and height of 20 cm. A layer of geotextile was pasted at the bottom of the pollution source chamber drilled with holes, in order to prevent the sand from entering into the pollution source chamber. A peristaltic pump was used to pump the contaminant with constant concentration into the pollution source chamber at the constant speed. Ten pressure taps were arranged at the bottom of the sand box, and each tap was connected with a transparent polyurethane pipe to measure the variation of the groundwater level. A 0.2 mm filter screen was set at the pressure tap, in order to prevent the sand from entering the pressure pipe and blocking the pipe. At the back of the sand tank, there were nine sampling ports in each layer. The distance between each two sampling ports was 20 cm. During the experiment, the sample was taken from the sampling port using a syringe. Then, the water sample was diluted and the electrical conductivity of the sample was measured by the DDS-307 conductivity meter. The concentration of Cl − can be calculated, because there is a linear relationship between the concentration of NaCl and electrical conductivity. Two sets of tests were considered for sea level rise. In the first experiment (Case 1), the left inland head and right sea level were fixed at 45 and 46 cm, respectively. In the second experiment (Case 2), the constant inland head and sea water level were maintained at 45 and 47 cm, respectively. The water level at different positions was observed at the measuring pressure holes during the experiment. The deionized water was used to make freshwater, and the density of freshwater was 1.0 × 10 3 kg/m 3 . The concentration and density of saltwater were 35.0 g/L and 1.02 × 10 3 kg/m 3 , respectively. The bright blue was added into the right seawater chamber, in order to show the seawater intrusion process. The saltwater with the concentration of 35.0 g/L was injected into the pollution source chamber. The carmine tracer with the concentration of 1.0 g/L was added into the saltwater to trace the movement of contaminant.
Each layer of the aquifer was filled with freshwater at the beginning of experiment. Then, the seawater flowed into the seawater chamber from the reservoir, and overflow was discharged through the overflow of the seawater chamber. Subsequently, the seawater gradually entered the water tank from the seawater chamber. The steady flow was formed on both sides. When the length of seawater wedge did not change in 5 min, the interface between seawater and freshwater reached steady state. Once the steady state was established, the peristaltic pump was open and the contaminant was injected into the pollution source chamber. The flow was set to 0.44 m/h, by controlling the peristaltic pump. At this time, the contaminant moved, influenced by the hydraulic gradient and density difference between seawater and freshwater. Finally, the experiment was finished when the contaminant plume intersected with the interface of seawater and freshwater in the coarse sand layer. The experimental process of Case 1 was similar to that of Case 2. During the experiment, the Canon digital camera (IXUS 175) was used to record the process of seawater wedge and contaminant migration. The salinity at each sampling port was monitored every 5 min. At the end of the experiment, the contaminant plume and saltwater wedge were analyzed by Matlab software.

Numerical Model and Procedure
Based on the finite difference model MODFLOW, the influence of density on groundwater flow was considered. A numerical simulation model of solute transport in groundwater flow was established by using SEAWAT [50,51]. The softwater was widely used in seawater intrusion and groundwater discharge e.g., [2,28,[52][53][54]. The coupling controlling equations of groundwater flow and solute transport with variable density were as follows: where ρ is the density of saline groundwater at a point in aquifer ; and R n is reaction term of chemical substance. The numerical simulation area was a heterogeneous and isotropic unconfined aquifer. During the setup of numerical simulation, the laboratory settings and parameters were followed as closely as possible. A two-dimensional vertical cross section was built, with the length of 220 cm and height of 55 cm. The domain of model was discretized into 24 layers and 5000 columns. A free surface boundary was defined on the upper part. A no-flow boundary condition was set on the bottom of the numerical model, because the bottom of the aquifer was bedrock in reality. For Case 1, the inland head and seawater level were set to 45 and 46 cm on the left and right boundaries, respectively. The left inland head and right seawater level were fixed at 45 and 47 cm for Case 2. The concentrations of saltwater and freshwater were set to 35 and 0 g/L, respectively. It was assumed that the contamination on the water table of the model was uniformly distributed at the distance between 1.04 and 1.08 m. The constant concentration of contamination was set to 35 g/L for Case 1 and Case 2. The injection time was at 55 and 50 min for Case 1 and Case 2, respectively. The simulated period was 110 min for Case 1 and 100 min for Case 2. The time step was set to 60 s for both cases.
In order to fit the observed head and salinity, the parameters of aquifers were calibrated repeatedly by a trial-and-error method. 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 the mean error between the calculated salinity values at the sampling ports was less than 10% of the maximum observed value variations. The horizontal hydraulic conductivity was assumed to be the same to its vertical value of the aquifer system. The estimated hydraulic conductivities were 4.6 × 10 −5 m/s, 2.9 × 10 −4 m/s, and 6.9 × 10 −4 m/s, respectively, for the fine sand layer, medium sand layer, and coarse sand layer, respectively. The other parameters used in the simulations were listed in Table 1.

Numerical Validation
The numerical simulations of variable density groundwater flow and the solute transport model reproduced the observed water, Cl − concentration, and contamination concentration in the transect. Here, the focus was on the Cl − concentration. Figure 2 shows the experimental versus the simulated Cl − concentrations variation with time at the four sampling points for Case 1 and Case 2. From the figure, one can see that the simulated and observed Cl − exhibited good correlation generally. The calculated correlation coefficient (R 2 ) ranged from 0.96 to 0.98, which indicated that the numerical model can successfully simulate the contaminant movement in the aquifer. From Figure 2, it can be seen that the observed and simulated values of Cl − concentration increased sharply at the initial time, and then they increased slowly and tended to be stable.
Here, the focus was on the Cl concentration. Figure 2 shows the experimental versus the simulated Cl concentrations variation with time at the four sampling points for Case 1 and Case 2. From the figure, one can see that the simulated and observed Cl exhibited good correlation generally. The calculated correlation coefficient (R 2 ) ranged from 0.96 to 0.98, which indicated that the numerical model can successfully simulate the contaminant movement in the aquifer. From Figure 2, it can be seen that the observed and simulated values of Cl concentration increased sharply at the initial time, and then they increased slowly and tended to be stable.
The Cl concentration values at the sampling ports No. 4″ and No. 5″ in the upper layer were higher than those at No. 4′ and No. 5′ in the middle layer. It indicated that the concentration of contaminant decreased gradually when the contaminant plume migrated downward in the aquifer. The migration of plume in the medium sand layer lagged behind that in the fine sand layer. From Figure 2a, one can see that the maximum value of Cl concentration at left No. 4″ was close to 35 g/L, which was higher than that at right No. 5″ in the upper layer. The value of Cl concentration at left No. 4′ was higher than that at right No. 5′ in the middle layer, which was similar to that in the upper layer. It was evident that the distribution of contaminant was affected by the flow from the sea to inland. Compared with Figure 2a, Figure 2b showed that the values of concentration at the sampling points (No. 4″, No. 5″, No. 4′, and No. 5′) in Case 2 were little smaller than those in Case 1. It demonstrated that the contaminant was diluted due to the increase of water flow velocity influenced by the sea level rising.

Seawater Intrusion Effected by Seawater Level Rise
In order to investigate the groundwater flow and salt transport under the influence of seawater level rising, simulations of seawater intrusion processes for different cases were conducted. Figure 3 reported the variation of observed data and simulated concentration color maps of the seawater wedge with time for Case 1. The comparison of Figure 3a,b reflected that the simulated toe lengths of the seawater wedge agreed with the observed ones well. The observed and simulated seawater wedge migrated inland with time, under the effects of the hydraulic gradient toward the inland and the density difference between saltwater and freshwater. The shape of the seawater wedge was parabolic. At the intial time, the migration velocity of the saltwater wedge was fast in the aquifer. The area increased in each layer of the aquifer. As for the hydraulic conductivity, the area and migration velocity of the seawater wedge increased from the top layer to bottom layer. Then, the migration velocity of the seawater wedge slowed down gradually, and it was close to zero when the time was 55 min. The interface of saltwater and freshwater reached equilibrium state, when the saltwater wedge did not move in the next five min. The toe position of 50% isoline of the seawater wedge represented the length of seawater wedge, which was wildly applied in depicting the saltwater wedge in the coastal aquifer, in which the dispersivity was small. As the invasion of seawater reaches the equilibrium status at the time of 55 min, the horizontal length of the 50% isoline of the observed seawater wedge can reach 142 cm. When the intruded seawater wedge was stable, the simulated length of the 50% isoline of seawater wedge was 144 cm, which was consistent with that of the observed seawater wedge.

Seawater Intrusion Effected by Seawater Level Rise
In order to investigate the groundwater flow and salt transport under the influence of seawater level rising, simulations of seawater intrusion processes for different cases were conducted. Figure 3 reported the variation of observed data and simulated concentration color maps of the seawater wedge with time for Case 1. The comparison of Figure 3a,b reflected that the simulated toe lengths of the seawater wedge agreed with the observed ones well. The observed and simulated seawater wedge migrated inland with time, under the effects of the hydraulic gradient toward the inland and the density difference between saltwater and freshwater. The shape of the seawater wedge was parabolic. At the intial time, the migration velocity of the saltwater wedge was fast in the aquifer. The area increased in each layer of the aquifer. As for the hydraulic conductivity, the area and migration velocity of the seawater wedge increased from the top layer to bottom layer. Then, the migration velocity of the seawater wedge slowed down gradually, and it was close to zero when the time was 55 min. The interface of saltwater and freshwater reached equilibrium state, when the saltwater wedge did not move in the next five min. The toe position of 50% isoline of the seawater wedge represented the length of seawater wedge, which was wildly applied in depicting the saltwater wedge in the coastal aquifer, in which the dispersivity was small. As the invasion of seawater reaches the equilibrium status at the time of 55 min, the horizontal length of the 50% isoline of the observed seawater wedge can reach 142 cm. When the intruded seawater wedge was stable, the simulated length of the 50% isoline of seawater wedge was 144 cm, which was consistent with that of the observed seawater wedge.  Figure 4 showed the variation of observed data and simulated results of the toe position of seawater wedge at various time for Case 2. The seawater wedge encroached into the inland with time, which was influenced by the density difference between saltwater and freshwater and the hydraulic gradient slope toward the sea. Similar to Case 1, the shape of the seawater wedge was parabolic. The migration velocity of seawater intrusion was fast initially, and then it slowed down gradually in Case 2. The migration velocity of the seawater wedge was close to zero at the time of 50 min, and the steady state of seawater intrusion was formed at this time in the three layers. The observed length of the 50% isoline of the seawater wedge was 167 cm, which was larger than that of Case 1. It was because the head at the seaside boundary and the flow velocity for Case 2 were higher than those for Case 1. It reflected that the rising sea level had a positively impact on the seawater intrusion process in the layered aquifer, which was consistent with the observation in Guo et al. [18]. At the time of 50 min, the simulated length of 50% isoline of the seawater wedge was 167 cm (Figure 4b), which was consistent with the observed value (Figure 4a). In general, the numerical simulation model can describe the variation of the seawater wedge well.  Figure 4 showed the variation of observed data and simulated results of the toe position of seawater wedge at various time for Case 2. The seawater wedge encroached into the inland with time, which was influenced by the density difference between saltwater and freshwater and the hydraulic gradient slope toward the sea. Similar to Case 1, the shape of the seawater wedge was parabolic. The migration velocity of seawater intrusion was fast initially, and then it slowed down gradually in Case 2. The migration velocity of the seawater wedge was close to zero at the time of 50 min, and the steady state of seawater intrusion was formed at this time in the three layers. The observed length of the 50% isoline of the seawater wedge was 167 cm, which was larger than that of Case 1. It was because the head at the seaside boundary and the flow velocity for Case 2 were higher than those for Case 1. It reflected that the rising sea level had a positively impact on the seawater intrusion process in the layered aquifer, which was consistent with the observation in Guo et al. [18]. At the time of 50 min, the simulated length of 50% isoline of the seawater wedge was 167 cm (Figure 4b), which was consistent with the observed value (Figure 4a). In general, the numerical simulation model can describe the variation of the seawater wedge well.

Contamination Transport Effected by Seawater Level Rise
The contaminate transport was simulated for both Case 1 and Case 2, when the seawater wedge was stable. Figure 5 showed the comparison between the observed and simulated contamination concentration at different times in the physical and numerical model for Case 1. The results showed that the numerical model could reproduce relatively well the overall shape of the plume, albeit the finger shape in the front of the plume observed in the laboratory was not captured in the simulation. The observed and simulated plume diffused mainly in the lateral and downward directions in the fine sand layer during the period of 55 and 70 min, and 70 and 85 min, respectively. The area of the contaminant plume increased linearly, and the shape of the plume was semi-elliptical. The maximum area of the contaminant plume was 0.02 m 2 , and the maximum depth of the contaminant plume was 0.13 m. At the time of 85 min, the contaminant entered into the medium sand layer and moved downward. The pollution plume became concave at the interface between the upper layer and middle layer, because the velocity of the contaminant plume migration increased gradually. The contaminant plume migrated in the medium sand layer at the time between 85 and 100 min. The variation rate of area in the medium sand layer was twice as much as that of the fine sand layer. The maximum area of the contaminant plume was 0.067 m 2 , and the maximum depth of the plume was 0.29 m. During the time of 100 and 110 min, the contaminant plume moved downward and inland in the coarse sand layer. The plume became concave at the interface between the middle layer and lower layer, due to the increasing of velocity of contaminant plume migration. The front of the pollutant flowed in the form of dominant flow, due to the heterogeneity in the coarse sand layer. The variation rate of the area in the coarse sand layer was 1.5 times of that in the middle sand layer, and the maximum plume area was 0.11 m 2 at the time of 110 min.

Contamination Transport Effected by Seawater Level Rise
The contaminate transport was simulated for both Case 1 and Case 2, when the seawater wedge was stable. Figure 5 showed the comparison between the observed and simulated contamination concentration at different times in the physical and numerical model for Case 1. The results showed that the numerical model could reproduce relatively well the overall shape of the plume, albeit the finger shape in the front of the plume observed in the laboratory was not captured in the simulation. The observed and simulated plume diffused mainly in the lateral and downward directions in the fine sand layer during the period of 55 and 70 min, and 70 and 85 min, respectively. The area of the contaminant plume increased linearly, and the shape of the plume was semi-elliptical. The maximum area of the contaminant plume was 0.02 m 2 , and the maximum depth of the contaminant plume was 0.13 m. At the time of 85 min, the contaminant entered into the medium sand layer and moved downward. The pollution plume became concave at the interface between the upper layer and middle layer, because the velocity of the contaminant plume migration increased gradually. The contaminant plume migrated in the medium sand layer at the time between 85 and 100 min. The variation rate of area in the medium sand layer was twice as much as that of the fine sand layer. The maximum area of the contaminant plume was 0.067 m 2 , and the maximum depth of the plume was 0.29 m. During the time of 100 and 110 min, the contaminant plume moved downward and inland in the coarse sand layer. The plume became concave at the interface between the middle layer and lower layer, due to the increasing of velocity of contaminant plume migration. The front of the pollutant flowed in the form of dominant flow, due to the heterogeneity in the coarse sand layer. The variation rate of the area in the coarse sand layer was 1.5 times of that in the middle sand layer, and the maximum plume area was 0.11 m 2 at the time of 110 min. Figure 6 showed the comparison between the transient experimental and numerical contaminant transport for Case 2. In all, the results showed that the transient movement of the plume was reasonably well simulated by the numerical model, although the finger shape in the front of the plume observed in the laboratory was not captured in the simulation. During the period of 50 and 70 min, the observed and simulated plume diffused laterally in the fine sand layer. The area of the contaminant plume increased linearly, and the shape of the plume was semi-elliptical, which was similar to Case 1. The maximum area and depth of the contaminant plume were 0.025 m 2 and 0.13 m, respectively. At the time of 70 min, the contaminant entered into the medium sand layer and moved downward. The pollution plume became concave at the interface between the upper layer and middle layer, because the velocity of the contaminant plume migration increased gradually. The contaminant plume migrated in the medium sand layer at the time between 70 and 90 min. The maximum area of the contaminant plume was 0.069 m 2 , and the maximum depth of plume was 0.29 m. During the time of 90 and 100 min, the contaminant plume moved downward and inland in the coarse sand layer. The plume became concave at the interface between the middle layer and lower layer, due to the increasing of velocity of contaminant plume migration. At the time of 100 min, the observed front of the pollutant was lower than that of the simulated one. Due to the heterogeneity in the coarse sand layer, the front of the pollutant flowed in the form of dominant flow. The variation rate of the area in the coarse sand layer was 2.0 times of that in the middle sand layer, and the maximum plume area was 0.13 m 2 at the time of 100 min. The end of the front of the contaminant plume migrated towards inland. Compared with the observed one, the numerical plume has a deviation toward the inland, because of the effect of seawater intrusion. Compared with Case 1, the area of the plume in Case 2 was larger than that in Case 1 in the same layer. It indicated that the velocity of diffusion of the contaminant plume increased as the sea level increased. As the velocity of diffusion increases, the time for the observed pollutant migration to the intersection between seawater and freshwater in Case 2 was shorter than that in Case 1.   Figure 6 showed the comparison between the transient experimental and numerical contaminant transport for Case 2. In all, the results showed that the transient movement of the plume was reasonably well simulated by the numerical model, although the finger shape in the front of the plume observed in the laboratory was not captured in the simulation. During the period of 50 and 70 min, the observed and simulated plume diffused laterally in the fine sand layer. The area of the contaminant plume increased linearly, and the shape of the plume was semi-elliptical, which was similar to Case 1. The maximum area and depth of the contaminant plume were 0.025 m 2 and 0.13 m, respectively. At the time of 70 min, the contaminant entered into the medium sand layer and moved downward. The pollution plume became concave at the interface between the upper layer and middle layer, because the velocity of the contaminant plume migration increased gradually. The contaminant plume migrated in the medium sand layer at the time between 70 and 90 min. The maximum area of the contaminant plume was 0.069 m 2 , and the maximum depth of plume was 0.29 m. During the time of 90 and 100 min, the contaminant plume moved downward and inland in the coarse sand layer. The plume became concave at the interface between the middle layer and lower layer, due to the increasing of velocity of contaminant plume migration. At the time of 100 min, the observed front of the pollutant was lower than that of the simulated one. Due to the heterogeneity in the coarse sand layer, the front of the pollutant flowed in the form of dominant flow. The variation rate of the area in the coarse sand layer was 2.0 times of that in the middle sand layer, and the maximum plume area was 0.13 m 2 at the time of 100 min. The end of the front of the contaminant plume migrated towards inland. Compared with the observed one, the numerical plume has a deviation toward the inland, because of the effect of seawater intrusion. Compared with Case 1, the area of the plume in Case 2 was larger than that in Case 1 in the same layer. It indicated that the velocity of diffusion of the contaminant plume increased as the sea level increased. As the velocity of diffusion increases, the

Sensitivity Analysis
The hydraulic conductivity, dispersivity, contaminate concentration, and heterogeneity were important parameters for solute transport in the coastal aquifer. Therefore, the sensitivities of

Sensitivity Analysis
The hydraulic conductivity, dispersivity, contaminate concentration, and heterogeneity were important parameters for solute transport in the coastal aquifer. Therefore, the sensitivities of saltwater contamination transported in the aquifer to the model parameters were analyzed, taking the Case 2 as an example. The area and velocity of the plume diffusion were discussed by changing the above parameters. It should be noted that the parameter values for analyzing the sensitivity were increased or decreased, whereas the values of the other parameters were not changed, which was shown in Table 1 Figure 6b, Figure 7a showed that the match to the observed contours of the saltwater contamination in the new model was worse, compared with the simulated ones in the Case 2. It was found that the larger length of the seawater wedge than in reality moved into the aquifer of the new model. The region of the transition zone in the front edge of the plume in the new model became larger, compared with that in the basic model. This was because the velocities were larger with higher hydraulic conductivities. The plume migrated more downward and landward, especially in the coarse sand layer. Increasing the hydraulic conductivity increased remarkably the area of the plume in the heterogeneous aquifer. The maximum area of the pollution plume simulated by the new model reached 0.143 m 2 , which was larger than that simulated in Case 2.
Sustainability 2020, 12, x FOR PEER REVIEW 11 of 17 was smaller than that of Case 2. Therefore, one can conclude that increasing or decreasing the hydraulic conductivities by the new models were not able to reproduce the spatial distributions of saltwater contamination concentration.

Effect of Dispersion
The sensitivity of the plume transport to the dispersivity was explored for two new models, due to the importance of dispersivity for solute transport. The variation of contours of the saltwater contamination was used for analyzing sensitivity with respect to the dispersivity. When the value of L α increased to 0.5 m and the ratio L T α α was 0.1, Figure 8a showed that more saltwater intruded into the aquifer of the new model than that of the basic case. The length of seawater in the new model was 0.13 m longer than that in Case 2. The simulated maximum vertical pollution depth was larger, compared with that in the basic model, which could be because high dispersivity induced widening of the transition zone. The area of the pollution plume simulated by the new model was 0.163 m 2 , which was larger than that of Case 2. It indicated that higher dispersion has great influence on plume transport. However, there was no change in the direction of pollutant migration.   Figure 7b showed that the match to the observed contours of the saltwater contamination was considerably worse, compared with the simulated ones in Figure 6b of Case 2. One can see that the smaller length of the seawater wedge than in reality migrated into the aquifer of the new model. Obviously, the maximum vertical pollution depth decreased when the hydraulic conductivities of the heterogeneous aquifer decreased, which indicated that the hydraulic conductivities for sensitivity analysis were too low. Additionally, with the lower hydraulic conductivities in each layer, the velocity of plume migration became slowly. However, the migration direction of the pollution plume remained downward and landward. The area of the pollution plume simulated by the new model increased slowly, and the maximum value of it reached 0.083 m 2 , which was smaller than that of Case 2. Therefore, one can conclude that increasing or decreasing the hydraulic conductivities by the new models were not able to reproduce the spatial distributions of saltwater contamination concentration.

Effect of Dispersion
The sensitivity of the plume transport to the dispersivity was explored for two new models, due to the importance of dispersivity for solute transport. The variation of contours of the saltwater contamination was used for analyzing sensitivity with respect to the dispersivity. When the value of α L increased to 0.5 m and the ratio α T /α L was 0.1, Figure 8a showed that more saltwater intruded into the aquifer of the new model than that of the basic case. The length of seawater in the new model was 0.13 m longer than that in Case 2. The simulated maximum vertical pollution depth was larger, compared with that in the basic model, which could be because high dispersivity induced widening of the transition zone. The area of the pollution plume simulated by the new model was 0.163 m 2 , which was larger than that of Case 2. It indicated that higher dispersion has great influence on plume transport. However, there was no change in the direction of pollutant migration.

Effect of Contaminant Concentration
The saltwater concentration was reduced from 35 to 20 g/L in the new model. Figure 9 shows that the fitting to the observed plume was worse, compared with that of the basic case in Figure 6b. One can see that the area of the pollutant plume and maximum vertical pollution depth decreased sharply, whereas the movement of the plume decreased slightly in the horizontal direction of the new model, compared with the simulated ones in Case 2. It was because the dispersion was too small. At the time of 100 min, the simulated area of the plume was 0.063 m 2 , which was smaller than that of the basic case. However, the direction of the pollutant plume migration was downstream and landward in the aquifer. Therefore, the range of the plume was sensitive to the concentration of contamination. When the value of α L decreased to 0.01 m and the ratio α T /α L was 0.1, Figure 8b showed that the area simulated by the new model could not match the observed ones at different times, compared with Figure 6b. The length of seawater in the new model was 0.093 m smaller than that in Case 2. It depicted that the area of pollution plume simulated by the new model was 0.086 m 2 , which was much smaller than that of Case 2. It indicated that decreasing the dispersion would lead to the slow migration velocity of plume. Overall, changes in dispersivity have effects on plume transport.

Effect of Contaminant Concentration
The saltwater concentration was reduced from 35 to 20 g/L in the new model. Figure 9 shows that the fitting to the observed plume was worse, compared with that of the basic case in Figure 6b. One can see that the area of the pollutant plume and maximum vertical pollution depth decreased sharply, whereas the movement of the plume decreased slightly in the horizontal direction of the new model, compared with the simulated ones in Case 2. It was because the dispersion was too small. At the time of 100 min, the simulated area of the plume was 0.063 m 2 , which was smaller than that of the basic case. However, the direction of the pollutant plume migration was downstream and landward in the aquifer. Therefore, the range of the plume was sensitive to the concentration of contamination.
sharply, whereas the movement of the plume decreased slightly in the horizontal direction of the new model, compared with the simulated ones in Case 2. It was because the dispersion was too small. At the time of 100 min, the simulated area of the plume was 0.063 m 2 , which was smaller than that of the basic case. However, the direction of the pollutant plume migration was downstream and landward in the aquifer. Therefore, the range of the plume was sensitive to the concentration of contamination.

Effect of Homogeneous Aquifer
In order to discuss the influence of heterogeneity on the contamination transport, numerical contamination simulation in the homogeneous aquifer was considered. The aquifer was assumed to be filled with fine sand based on Case 2. The value of K of the whole aquifer was taken as 4.0 m/d, and the other parameters were fixed as shown in Table 1. Figure 10 showed the simulated contamination concentration distributions of the cross section in the homogeneous aquifer for the new model. The steady state length of the seawater in the new model was smaller than that in Case 2. At the time of 70 min, the shape of the pollution plume was semi-elliptic, which was the same to that of the basic case. Then, the plume migrated downward and laterally during the period between 70 to 100 min, compared to that of the basic case. The vertical migration of the plume decreased, nevertheless, the movement in the horizontal direction increased slightly for the new model. This was because the vertical velocity of the new model was smaller than that of the basic case. There was no concave observed in the plume during the process of pollutant transport, compared to Figure 6b. At the time of 100 min, the maximum area and vertical pollution depth of the pollutant plume reached 0.11 m 2 and 0.25 m in the new model, which were smaller than those simulated by the basic case. Therefore, one can conclude that the infiltration depth and range of the contaminant plume in the heterogeneous aquifer were greater than those in the homogeneous aquifer of the actual beach.
was because the vertical velocity of the new model was smaller than that of the basic case. There was no concave observed in the plume during the process of pollutant transport, compared to Figure 6b. At the time of 100 min, the maximum area and vertical pollution depth of the pollutant plume reached 0.11 m 2 and 0.25 m in the new model, which were smaller than those simulated by the basic case. Therefore, one can conclude that the infiltration depth and range of the contaminant plume in the heterogeneous aquifer were greater than those in the homogeneous aquifer of the actual beach.

Conclusions
This paper presented for the first time analyzing the influence of sea level rise on contaminate transport in the heterogeneous unconfined aquifer of the coastal zone. Two sets of laboratory experiments were tested considering the rising sea level (the difference between sea level and inland head was 1 cm for Case 1 and 2 cm for Case 2). The variation of seawater intrusion and contaminant transport with time were observed during both sets of experiments. The numerical model SEAWAT

Conclusions
This paper presented for the first time analyzing the influence of sea level rise on contaminate transport in the heterogeneous unconfined aquifer of the coastal zone. Two sets of laboratory experiments were tested considering the rising sea level (the difference between sea level and inland head was 1 cm for Case 1 and 2 cm for Case 2). The variation of seawater intrusion and contaminant transport with time were observed during both sets of experiments. The numerical model SEAWAT was used for validation. The numerical results of the model provided matched well with the experimental data for both the toe length of seawater wedge and the shape of the contaminant plume for the rising sea level in the heterogeneous unconfined aquifer.
The observed and simulated values of Cl − concentration at the sampling points increased sharply at the initial time, and then they increased slowly and tended to be stable. The Cl − concentration in the medium sand layer was smaller than that in the fine sand layer, which indicated that the migration of plume in the middle layer lagged behind that in the upper layer. The values of concentration at the sampling points in Case 2 were a little smaller than those in Case 1, because the contaminant was diluted due to the increase of water flow velocity influenced by the rising sea level. The observed and simulated seawater wedge migrated inland with time, under the effects of the hydraulic gradient toward the inland and the density difference between saltwater and freshwater. The steady state length of the 50% isoline of the seawater wedge was 167 cm in Case 2, which was larger than that of Case 1. It was because the head at the seaside boundary and the flow velocity for Case 2 were higher than those for Case 1.
In both cases, the observed and simulated plume diffused laterally and downward in the fine sand layer and medium sand layer. The contaminant plume moved downward and inland in the coarse sand layer. The front of the pollutant flowed in the form of the dominant flow, due to the heterogeneity of the aquifer. At the interface between each two layers, the plume became concave due to the increasing velocity of contaminant plume migration. The maximum area of plume in Case 2 was 0.13 m 2 , larger than that in Case 1, which indicated that the velocity of diffusion of the contaminant plume increased as the sea level increased. As the velocity of diffusion increased, the time for pollutant