Abstract

Laboratory experiments were conducted in a vertical, two-dimensional, rectangular flow tank, simulating the response of a phreatic coastal aquifer to a sea tide. Imposed sinusoidal fluctuations of the saltwater level at one side of the flow tank caused three types of fluctuations: (a) hydraulic head throughout the aquifer, (b) saturation degree within the capillary fringe, and (c) salt concentration surrounding the freshwater-saltwater interface (FSI), all recorded by head, saturation, and salinity sensors, respectively. Significant time lags were observed both in the saturation degree within the unsaturated zone and in the salinity within the FSI. All measured values, recorded by the three types of sensors, were simulated and reproduced using a numerical model. The calibrated model was used for mapping the time lags throughout the aquifer. It was found that the time lag of saturation fluctuations within the unsaturated zone increased upward from the groundwater level as the unsaturated hydraulic conductivity decreased. Similarly, the time lag of salinity fluctuations within the FSI increased downward, with distance from the groundwater level. We interpret the low hydraulic conductivity at the capillary zone as the source of attenuation of both saturation and salinity, because both are controlled by the vertical advection of the whole freshwater body. This advection is significantly slower compared to the dynamics of pressure diffusion. The uniqueness of this study is that it provides quantitative data on the attenuation at the capillary zone and its effect on the salinity time lag in coastal aquifer systems.

1. Introduction

Seawater intrusion into coastal aquifers is of a global concern due to salinization of production wells (e.g., [1, 2]). The transition zone between the fresh groundwater that flows seaward and the saline seawater beneath that penetrates inland is referred to as the freshwater-saltwater interface (FSI). The two water bodies are in a dynamic equilibrium that is influenced by hydrodynamic processes such as sea tide oscillations and pumping [3, 4].

Sea tide oscillations influence both the groundwater level (GWL) and the location of the FSI in coastal aquifers. They induce groundwater level fluctuations with a periodicity similar to tidal periodicity [5, 6], as well as salinity fluctuations within the FSI [4, 7, 8]. GWL fluctuations in an unconfined aquifer subjected to low-frequency tidal oscillations are significantly affected by the capillary effect. Parlange and Brutsaert [9] showed that the capillary pressure above the water table varies as a consequence of GWL fluctuations. Nielsen and Perrochet [10] and Werner and Lockington [11] found a significant effect of the capillary fringe and the hysteresis phenomenon, on GWL fluctuations in an unconfined aquifer subjected to tidal oscillations. Cartwright et al. [12] showed the influence of water table shallowness on water table wave dispersion. Kong et al. [13, 14] used analytical solutions and numerical models to show that unsaturated flow above the water table influences the amplitude of the GWL fluctuations and the lag of GWL fluctuations behind the sea level.

The depth and width of the FSI are directly influenced by the groundwater head fluctuations; thus, the FSI itself fluctuates due to sea level oscillations as well [1, 7, 15]. Ataie-Ashtiani et al. [3] showed theoretically that the sea tide widens the FSI and forces the sea water to intrude further inland. Kim et al. [8] found that the FSI and the water table in Kimje, Korea, fluctuate with a tidal periodicity. Levanon et al. [16] demonstrated that the GWL and the FSI fluctuations in the coastal aquifer of Israel are composed of the same tidal constituents as the tide wave in the Mediterranean Sea. Heiss and Michael [17] conducted a field study and numerical model and showed that the size of the mixing zone varied over the spring-neap tidal cycle. They also showed that tidal amplitude significantly affects the flow of fresh and saline groundwater toward the sea. Pool et al. [18] quantify tidal impacts on solute mixing and spreading in the FSI using a numerical model and showed that the key parameters controlling the shape and location of the FSI are the tidal amplitude, the tidal period, and the hydraulic diffusivity.

Levanon et al. [19, 20] suggested that the tidal effect on the groundwater system comprises two main processes: (1) tidal forcing at the sea floor boundary and (2) attenuation at the groundwater level due to capillary effect. They showed that even though the pressure head wave propagation into the aquifer is relatively fast, the actual movement of the freshwater body, which is reflected by fluctuations of GWL and the salinity in the FSI, is slower since it is controlled by the unsaturated hydraulic conductivity within the capillary fringe. Up to date, this thesis was shown to occur by a field study, in which the horizontal time lags of both head and salinity were measured across the aquifer using wells. However, the actual delay within the unsaturated zone and the vertical time lag that progress downwards all the way to the FSI was not shown in the field. This can be examined in a controlled laboratory setup, where it is easier to quantify the vertical time lags and unsaturated processes.

The use of physical models provides an ideal setup for conceptualization, visualization, and measurement of groundwater flow and solute transport through an aquifer [21]. Physical models have been used in a variety of groundwater studies for describing different phenomena, such as density-driven groundwater flow, saltwater intrusion, and transport of solutes and contaminants (e.g., [2225]). Physical models were also used to show the effect of sea tides on the residence time of a plume that originated at the sea that infiltrates into the aquifer downwards during high tides and back to the sea during low tides [2628]. However, in these experiments, the positions of the FSI were not established and monitored. The physical model in this study provides a better understanding of the dynamics in phreatic coastal aquifers subjected to sea tide influence including the FSI, water level, and saturation.

The objective of this paper is to investigate the importance of the unsaturated flow properties within the capillary fringe on the tidal dynamics of both the FSI and the water table fluctuations by controlled laboratory experiments and numerical modeling. While in previous studies [19, 20] only the horizontal time lag induced by the sea tide was analyzed at the field scale, in this study, the vertical development of the time lag is displayed; thus, a more complete picture of the tidal dynamics in coastal aquifers is achieved.

2. Methods

2.1. Laboratory Setup

Laboratory experiments were conducted in a vertical, rectangular (100 cm long, 50 cm high, and 3 cm wide) flow tank made of Plexiglas, simulating a phreatic cross-shore coastal aquifer (Figure 1). The narrow width of the flow tank ensures that the problem examined could be defined as two-dimensional, with negligible along-shore processes. It is divided into three distinct chambers: a central flow chamber that serves as an aquifer and contains the porous medium and two side chambers filled with water, which define the boundary conditions during the experiments. The left chamber represents the saltwater boundary at the seaside, and the right chamber represents the inland freshwater boundary of the regional aquifer. These side chambers are separated from the main chamber by a fine net that prevents the passage of granular material.

The initial and boundary conditions during the experiment are defined by the water levels at the side chambers. At the right boundary, freshwater is pumped from the supply reservoir into an overflow bottle connected to the right chamber. At the left boundary, the saltwater is pumped from another supply reservoir directly into the left chamber. The water level at the right boundary is higher than that at the left one; thus, an initial hydraulic gradient is formed, and the water flows through the porous medium from right to left simulating a seaward flow. The overflow at the left boundary is diluted sea water as a result of mixing between the freshwater and the saltwater which occurs along the flow paths in the FSI. The overflow pipe at the left chamber is connected to a computer-controlled engine that enables quasitidal fluctuations of the saltwater at the left boundary, which define the boundary condition at the seaside. The amplitude and the wavelength of the sinusoidal tidal wave are definable by the engine (changing vertical velocities and the time intervals for each velocity). The tidal wave used in the experiments is harmonic with constant amplitude. During all the experiments, the water levels at the boundaries were kept lower than 35 cm (Figure 1), in order to enable an unsaturated zone above the water table and to avoid problems that truncation of this area may create [29]. The freshwater in the experiment is colorless tap water, and the saline water is prepared by dissolving NaCl in tap water and adding commercial red food dye (AmeriColor Ltd.). The density of the saline water was 1.027 gr cm-3, the same as the average density of the eastern Mediterranean Sea. The sand used in the experiments is silica grains which were sieved to a diameter range of 800-1,300 μm. In preexperiment treatment, it was washed with distilled water to remove dust and clay minerals and with HCl acid to remove the oxide coatings of the quartz grains. The horizontal permeability was calculated from a simple steady state experiment where the hydraulic heads were constant at both ends and a freshwater was flowing from the high hydraulic head at the left boundary to the low hydraulic head at the right boundary, using the equation . The porosity was measured by dividing the weight of the sand to the volume of the tank. The values of the parameters are given in Table 1.

The monitoring setup during the experiments includes several measurement systems (Figure 1(b)): (1) 168 electrodes placed at the back side of the flow tank for in situ voltage measurements, which are equivalent to the salinity of the water. The data from the electrodes is collected using LabVIEW software (National Instruments Ltd.); (2) two water level sensors (PA-36XW, Keller Ltd.) located at the side chambers; (3) four water content sensors (EC-5; Decagon Devices Inc.) placed in the capillary zone above the water level; and (4) six piezometers located a few centimeters below the water table at different distances from the left boundary. In four of these piezometers (P1-P4 in Figure 1), a thin film of red oil is added on top of the water column in order to allow easier monitoring of the changes in the water level. A high-resolution video camera (50 frames per second) records the piezometers, and the data is converted to numerical data using image processing. In addition, a pressure sensor measured the barometric pressure in the laboratory during the experiment for barometric compensation of the pressure measurements at the boundaries.

2.2. Estimation of the Unsaturated Zone Characters

Unsaturated flow is an important factor in coastal aquifers under tidal influence (e.g., [11, 13, 30]). To evaluate the unsaturated properties of the sand used in the flow tank experiments, a preexperiment was conducted, in which the volumetric water content was measured under different conditions of capillary tension. Evaluation of the unsaturated characters of the sand is important for a few reasons. First, it helps plan the location of the water content sensors in the experiments. Second, according to the retention curve, the experiment is designed so that a fully unsaturated zone will develop, and third, it is important for the calibration of the numerical model.

Two sensors of water content were located just above the water table in a wide column full of sand. An overflow bottle placed on an elevator was connected to the bottom of the column, which enabled changes of its water level. The water level in the tank was lowered stepwise until residual water content was measured. Then, the water level was raised gradually, in an opposite process, to its initial height. In each step, the system was static for 10-15 minutes to allow full drainage or filling of the porous medium. The retention curve of the drying process is different from that of the wetting process due to hysteresis (Figure 2). The air entry value for the drying process is located at ~5 cm above the water table, and for the wetting, it is located at ~3 cm. The residual water content of the sand is relatively low, only 0.025, due to the high level of the sieving. The saturated water content (equal to the porosity) is 0.38 for the drying process. However, for the wetting process, it was found to be 0.35, probably as a consequence of some trapped air.

The retention curve is quantitatively described by [31]: where (-) is the volumetric water content, (L) is the capillary tension, (-) and (-) are the saturated and residual volumetric water content, respectively, and (L-1), (-), and (-) are empirical constants, where .

By fitting the retention curve (equation 1) to the experimental results, the empirical unsaturated parameters for the drying and wetting processes were calibrated. The values for are 11.8 m-1 and 18 m-1 and for are 9 and 7.5 for the drying and wetting processes, respectively (Figure 2). These values were used in the numerical model (Table 1).

2.3. Numerical Model

Numerical simulation of the laboratory experiments is conducted using the FEFLOW code, a finite-element simulator that solves the coupled variable density groundwater flow and solute transport equations [32]. In the case of tide-induced fluctuations, both saturated and unsaturated flows are involved. The general equation for saturated-unsaturated flow, derived from the conservation of fluid mass, for variable densities is [32, 33] where (-) is the saturation degree; (L-1) is the specific storage coefficient; (-) is the effective porosity; (L) is the pressure head; (L) is the fresh groundwater hydraulic head; (T) is the time; (L T-1) is the hydraulic conductivity; , where (-) is the relative viscosity of the fluid, (M L-1 T-1) is the viscosity of water, and (M L-1 T-1) is the water viscosity at the standard state; , where (-) is the relative density, (M L-3) is the fluid density, and (M L-3) is the density at the standard state; and is the vertical vector unit [0,0,1].

The conservation of solute mass may be written as [1] where (M L-3) is the fluid concentration, (L2 T-1) is the dispersion-diffusion coefficient, and (L T-1) is Darcy’s velocity of the fluid. The coupling between the solute transport equation and the flow equation is through equations of state for the fluid’s density and viscosity.

The dimensions of the numerical cross-section are the same as the laboratory flow tank—a length of 100 cm and a height of 50 cm. A mesh with 228,090 elements and 114,752 nodes is used with refinement near the FSI in order to increase the stability of the numerical scheme and the accuracy of the model (Figure 3). The hydraulic heads at the left and right boundaries are defined as a tidal time series, based on the water level fluctuations that were measured during the experiment at the side chambers. The equation which defines these time series is where is the hydraulic head (L), is time (T), (L) is the initial hydraulic head, (L) is the tidal amplitude, and (Rad T-1) is the angular frequency of the tidal wave. The amplitudes of the left and right boundaries are 3 cm and 0.4 cm, respectively, and the wavelength is 15 minutes, in accordance with the laboratory experiment. At this frequency, the hysteresis effect is very small [34]. In addition, no-flow boundaries are defined at the cross-section base. Salinity on the sea side is 39,000 mg L-1 (based on the density and the temperature of the saline water in the laboratory), and on the right boundary is 20 mg L-1. Table 1 summarizes the parameters used in the numerical model. The calibration process is detailed below.

2.4. Cross-Correlation

In order to quantitatively analyze the time lags between sea level fluctuations, salinity, water level, and saturation in the aquifer, a cross-correlation analysis was applied on both laboratory and numerical model results. Such analysis provides information about the correlation between input and output time series, along with its significance. Details about the application of the cross-correlation analysis to coastal aquifer systems can be found in Levanon et al. [20].

3. Results

3.1. Laboratory Experiments

Snapshots from one of the tidal cycles are presented in Figure 4, along with the comparable numerical results. The FSI moves forward and backward under the influence of the tidal fluctuations at the left boundary (see Supplementary Materials (available here) for a video of the full dynamic process of the experiment). At its maximum penetration, the FSI toe reaches 39 cm from the left boundary (Figure 4(b)), and at its minimum, it reaches 30 cm (Figure 4(d)). The maximum and minimum penetration distances of the FSI toe lag a few minutes behind the high and low tides (i.e., maximum and minimum water levels at the left boundary, Figures 4(a) and 4(c)). The significant difference between the salinity distributions at high and low tides is in the contact area with the left boundary, where at the high tide the saline water is higher than at the low tide.

3.1.1. Horizontal Effect

With sampling frequency of 20 msec, no time lag has been detected between the water level at P1, which is located 18 cm from the left boundary, and the water level at P4, which is located 80 cm from the left boundary (Figure 5). This is attributed to the fast pressure head wave propagation into the aquifer which is controlled by the high diffusivity () of the aquifer [20]. There is also no lag between the salinity at electrodes which are at the same depth but at different distances from the left boundary (electrodes #7 and #8 in Figure 5). However, the fluctuations of the salinity lag significantly after those of the water level. This lag between the water level and the salinity is related to vertical flow which is considered next.

3.1.2. Vertical Effect

On the other hand, there is a significant lag between sensors that are vertically aligned. Within the unsaturated zone, above the water table, the upper water content sensor (WC1) lags behind the lower water content sensor (WC2) and they both lag behind the sea level fluctuations (Figure 6). The salinity at the FSI (electrode #28) lags behind both water content sensors. All these sensors are aligned vertically 18 cm from the left boundary. This vertical time lag trend exists also within the FSI where the lower electrode #9 also lags behind the upper electrode #13 (Figure 7). The mechanism for the vertical lag in saturation and salinity will be further discussed later.

3.2. Numerical Model

The hydraulic parameters (Table 1) used in the numerical simulations were calibrated to achieve a reasonable fitting between the experimental and modeled salinities (Figures 4 and 8). Even in the controlled environment of this experiment, there is some heterogeneity in the packing of the sand and some uncertainties in the measurements. Therefore, there is no attempt here to numerically reproduce the laboratory experiments exactly. Nevertheless, it helps to understand the mechanism responsible for the lags of the saturation and salinity. Based on the numerical simulation results, cross-correlation analysis was used to evaluate the correlation and the time lag between the sea level as an input and salinity and saturation as outputs.

In locations where the salinity and saturation are influenced by the sea tide, the lag after the sea level fluctuations is calculated using cross-correlation (Figure 9). Changes in salinity are observed only close to the FSI, and changes in saturation are observed only above the water level. In these locations, it is possible to calculate time lags. Because of the specific Van Genuchten parameters in this experiment, the significant drop in the saturation occurs about 9 cm above the water level. As a result, the unsaturated hydraulic conductivity also drops in this elevation and causes a significant lag for water content fluctuations. This lag of flow is transmitted to deeper parts of the aquifer and causes a lag of the salinity fluctuations around the FSI (Figure 9). The horizontal lag of the system to the imposed horizontal boundary condition is negligible because the pressure signal is transmitted rapidly through the cross-section whereas the vertical lag of the saturation and salinity is significant because it is related to the actual advection of the water.

4. Discussion

Usually, the use of a small-scale laboratory setup has a limitation due to the downscaling problem in time and space. However, such a laboratory physical setup can provide actual results that can verify the outcome of hydrological simulation in a place where no field data can be obtained for such verification. Thus, the laboratory results can have great advantages which help in understanding the hydrodynamic processes involved in tidal influence on the coastal aquifer. As mentioned before, the effect on the vertical scale could only be obtained at the lab experiment and not in the field.

The results of this study confirm the mechanism hypothesized by Levanon et al. [19, 20] for sea tide influence on the coastal groundwater system, including both the FSI and the GWL. The tidal effect on the groundwater system comprises two main processes: tidal forcing at the sea floor boundary and attenuation at the groundwater level due to capillary effect. This mechanism was previously proven at the horizontal scale by field monitoring and numerical simulations. The horizontal pressure head wave propagation into the aquifer is fast; thus, it could be detected only at the field scale. Indeed, at the laboratory scale, there is no measured horizontal time lag between the boundary level fluctuations and the GWL (Figure 5). The current study verifies the hypothesis through controlled laboratory experiments at the vertical scale by the time lags at both the unsaturated and the saturated zones. The vertical tide-induced fluctuations of the groundwater level create an actual movement of the whole water body in a phreatic aquifer. Within the capillary zone, these fluctuations are attenuated by the relatively low hydraulic conductivity. The exact extent of the attenuation is dictated by the unsaturated hydraulic conductivity which is related to the retention curve of the matrix. The specific properties of sand used in this study (Figure 2) caused time lags of the saturation degree up to nine cm above the water level (Figure 9). At the same time, the salinity at the FSI is fluctuating with the fluctuation of the whole water body. This fluctuation is controlled by an advection wave (as opposed to the horizontal pressure wave), and therefore, the lag increases with depth (Figure 9).

This study only examined an aquifer with a vertical shoreline and no hysteresis. A sloping shoreline will likely complicate patterns and behaviors of groundwater flow and disrupt the sinusoidal tide due to the asymmetry of the tidal infiltration/draining process for a sloping beach [6, 17, 35]. It is much easier for the water to infiltrate into a sloping beach at high tide than to drain away at low tide. Hysteresis of wetting and draining water level has a similar effect, a strong asymmetry in the moisture content response to the symmetric forcing condition. For periods less than 15 minutes, the dynamics become symmetric and nonhysteretic. [34]. To avoid these complications in this study, we chose a vertical shoreline and a period of 15 minutes for the tidal wave.

5. Conclusions

Laboratory experiments conducted in a vertical, two-dimensional, rectangular flow tank exhibit a phreatic coastal aquifer. Sea tide effects on both the FSI and the GWL are analyzed in the laboratory along with numerical models that simulate these experiments. Significant time lags are observed both in the saturation degree within the capillary zone and in the salinity within the FSI. The low hydraulic conductivity in the capillary zone attenuates the movement of the water table, which controls the movement of the freshwater aquifer and causes the saturation and salinity to lag behind the fluctuating sea. The time lag of the saturation is the highest at the top of the capillary zone where the saturation and permeability is the lowest. The time lag of the salinity within the FSI increases with depth with distance from the capillary zone.

Most previous field studies in coastal aquifers neglect the dynamic nature of the attenuation caused by the nonlinearity at the capillary zone. However, these dynamics may affect significantly the transport, mobility, and removal of chemicals and their subsequent discharge to coastal waters. The attenuation at the capillary zone and its effects on the entire coastal aquifer revealed here provide new insight into the complexity, intensity, and time scales of tidal effects in coastal aquifers.

Data Availability

All data are available on request at the following e-mail address: [email protected].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was financed by the Israeli Water Authority and the Israel Science Foundation, grant number 942/14. We want to thank I. Simchayoff and A. Podelko from the Hebrew University of Jerusalem who helped in planning and developing the laboratory system.

Supplementary Materials

A video of the laboratory experiment showing the full dynamic process of the experiment. (Supplementary Materials)