Transience of seawater intrusion and retreat in response to incremental water‐level variations

This paper provides for the first time an experimental study where the impact of sea‐level fluctuations and inland boundary head‐level variations on freshwater–saltwater interface toe motion and transition zone dynamics was quantitatively analysed under transient conditions. The experiments were conducted in a laboratory flow tank where various (inland and coastal) head changes were imposed to the system and the response of the key seawater intrusion parameters was analysed with high spatial and temporal resolution. Two homogeneous aquifer systems of different grain size were tested. The numerical code SEAWAT was used for the validation. The results show that in cases of sea‐level variations, the intruding wedge required up to twice longer time to reach a new steady‐state condition than the receding wedge, which thereby extend the theory of timescale asymmetry between saltwater intrusion and retreat processes in scenarios involving sea‐level fluctuations. The intruding and receding rates of the saltwater wedge were respectively similar in the scenario involving sea‐level and the freshwater‐level changes, despite change in transmissivity. The results show that, during the intrusion phase, the transition zone remains relatively insensitive, regardless of where the boundary head change occurs (i.e., freshwater drop or sea‐level rise) or its magnitude. By contrast, a substantial widening of the transition zone was observed during the receding phase, with almost similar amplitude in the scenario involving a rise of the freshwater level compared with that caused by a drop of the saltwater level, provided that an equivalent absolute head change magnitude was used. This transition zone widening (occurring during saltwater retreat) was greater and extended over longer period in the low hydraulic conductivity aquifer, for both freshwater‐level rise and sea‐level drop scenarios. The concentration maps revealed that the widening mechanism was also enhanced by the presence of some freshwater sliding and into the wedge during saltwater retreat, which was thereafter sucked upward towards the interface because of density difference effects.


| INTRODUCTION
For populations living in coastal zones, which represents about 70% of the world's population (Webb & Howard, 2011), groundwater represents the main, in some cases the only, source of water supply. In such areas, groundwater is particularly susceptible to degradation given its proximity with seawater as well as the intense water demand associated with ever-increasing population densities. The subsequent overabstraction of groundwater to satisfy an increasing demand ultimately leads to the landward incursion of oceanic saline water into fresh coastal groundwater commonly known as seawater intrusion (SWI). Because climate change is expected to aggravate in the near future and thus cause further detrimental effects of SWI, it has become of great importance to further advance the understanding of underlying processes related to SWI for effective water resource management . SWI is generally characterized by quantitative indicators, which include external metrics that delineate the outward boundary of the saline plume and the freshwater-saltwater interface. The external SWI metrics include the horizontal extent of intrusion, often called the toe length, and the vertical extent along the coastline boundary (Goswami & Clement, 2007;Sebben, Werner, & Graf, 2015). The freshwater-saltwater interface is a specific characteristic feature of coastal aquifers that occurs as a transition zone (TZ) where the salt concentration varies gradually from that of the freshwater to that of the saltwater in the seaward direction. These parameters outline the wedge-like shape of the plume that the intruding saltwater tend to form while encroaching homogeneous isotropic coastal aquifers, albeit aquifer heterogeneity may, to various degrees, severely affect this highly idealized shape.
In addition to unreasonable extraction of groundwater, climate change also greatly enhances the saltwater intrusion mechanism in coastal aquifers (Ferguson & Gleeson, 2012;Oude Essink, van Baaren, & de Louw, 2010;Sherif & Singh, 1999). The reduction in freshwater recharge, resulting from changes in rainfall patterns, and sea-level rises (SLRs) are the two climate change-induced factors that cause a disruption of the hydraulic gradients between the land and the sea, subsequently leading to more pronounced SWI. The reduction of freshwater recharge into an aquifer induces a direct decrease of the magnitude of the seaward groundwater discharge, which leads to a lesser resistance exerted on the intruding saline water. The rate of freshwater recharge into the aquifer and the magnitude of seaward groundwater discharge are considered the major controlling factors affecting SWI (Ketabchi, Mahmoodzadeh, Ataie-Ashtiani, & Simmons, 2016).
It is often difficult to derive a clear understanding of the mechanisms affecting SWI directly from field-based investigation . The challenge of measuring and quantifying coastal aquifer hydrodynamics and SWI in field sites has promoted the use of laboratory and numerical modelling tools to gain a valuable insight into SWI response to various geological and/or hydrological stresses, such as (a) change in seaward freshwater discharge resulting from fluctuations at the inland head boundary ( Robinson, Ahmed, & Hamill, 2016;Robinson, Hamill, & Ahmed, 2015) in head-controlled systems or from variations of the regional freshwater flux (Chang & Clement, 2012;Stoeckl & Houben, 2012;Stoeckl, Houben, & Dose, 2015) in flux-controlled systems and (b) SLR (Hussain & Javadi, 2016;Morgan, Bakker, & Werner, 2015;Morgan, Stoeckl, Werner, & Post, 2013). Goswami and Clement (2007) provided a thorough quantitative analysis of the steady state and transient toe response to freshwater head changes to provide an improved benchmark for numerical models. Chang and Clement (2012) examined the difference in the toe response to variations in freshwater flux resulting from regional and areal flow changes. Their study highlighted the difference of timescale associated with saltwater intrusion and retreat and provided explanation of the difference based on flow velocity field analysis.
Although these aforementioned studies have emphasized the impact of changes in freshwater flow on saltwater intrusion dynamics, they do not include a quantitative comparison of the toe migration rate in response to inland and coastal head changes. The timescale of toe response was the primary focus of Lu and Werner (2013), whereby the timescale of saltwater intrusion and retreat was related by simple linear empirical correlations to the boundary head difference and distance of toe response, respectively. Although their study examined the timescale of toe response to inland and coastal head changes, their analysis was based on a confined aquifer system, where the saturated thickness remains constant in response to SLR. This study aims to provide a quantitative analysis of the intruding and receding saltwater wedge migration rates in response to inland and coastal head fluctuations in unconfined aquifer systems.
The adverse impact of SLR on SWI has been well documented in the literature (Chang, Clement, Simpson, & Lee, 2011;Ketabchi et al., 2016;Ketabchi, Mahmoodzadeh, Ataie-Ashtiani, Werner, & Simmons, 2014;Michael, Russoniello, & Byron, 2013;Watson, Werner, & Simmons, 2010;Werner & Simmons, 2009). Hussain and Javadi (2016) showed that SLR significantly enhances the inland progression of saline water and caused a substantial reduction of the available fresh groundwater resources, particularly for flat-shoreline aquifer systems. The type of inland boundary condition was found to play a major role in determining the vulnerability of a coastal aquifer to SWI induced by SLR. In physical terms, substantially larger saltwater intrusion lengths have been observed in head-controlled aquifer systems compared with those observed in flux recharge-controlled systems (Michael et al., 2013;Rasmussen, Sonnenborg, Goncear, & Hinsby, 2013;Werner & Simmons, 2009). Several studies gave a special attention to the overshooting mechanism (Ketabchi et al., 2016;Morgan et al., 2013;Morgan et al., 2015;Watson et al., 2010), a process by which the saltwater intrusion length extends transiently beyond the final equilibrium position as a result of SLR, thereby questioning the steady state as being the actual worst case scenario.
Investigations on the impact of sea-level fluctuations on the dynamics (widening and narrowing) of freshwater-saltwater TZ are however very scarce and often limited to numerical modelling works. Lu, Kitanidis, and Luo (2009) explored TZ behaviour under transient flow condition using numerical simulations in a scaled-tank model and field model. They concluded that the widening process was primarily controlled by kinetic mass transfer combined with TZ motion induced by water-level fluctuations. Lu and Luo (2010) further demonstrated the sensitivity of TZ dynamics to the flow velocity.   Figure 1 shows a schematic diagram of the flow tank. The flow tank was subdivided into three distinct compartments: a central chamber and two circular reservoirs located at either side.
The central chamber was filled with clear glass beads from Whitehouse Scientific® to simulate the porous medium. Two fine mesh acrylic screens with aperture diameter of 0.5 mm were located at each side of the central chamber to separate the porous medium from the side reservoirs. The apertures of these meshes were small enough to contain the glass beads and sufficiently large to allow circulation of water flowing from the side reservoirs.
The two reservoirs were used to provide constant-head boundary condition at either side of the porous media. The left and right side reservoirs were filled with freshwater and saltwater, respectively.
The freshwater reservoir was filled with cold tap water. Adjustable overflow outlets were placed within each reservoir to maintain the water levels by draining any excess water to waste. Each of the overflow outlets was connected to a flexible hose that was plugged at the bottom of each reservoir as shown in Figure 1. The saltwater filling the right side reservoir was sourced from a 200-L saline water solution prepared prior the experiments. The saltwater solution was prepared by dissolving food grade into freshwater at a concentration of 36.16 g/L to achieve a density of 1,025 kg/m 3 .
The density of the saltwater solution was monitored using a hydrometer H-B Durac plain-form polycarbonate and through manual measurements using mass/volume ratio. In order to distinguish it from the freshwater, red food colour was added to the saltwater at concentration of 0.15 g/L. The scale used for the measurement of salt and dye quantities was accurate to 0.01 g. For each set of experiments, the saltwater was primarily sourced from the initial 200-L batch in order to reduce any possible risks of density or colour discrepancies between each experiment to a minimum. The water levels in both freshwater and saltwater reservoirs were monitored using ultrasonic sensors Microsonic mic+25/DIU/TC with ±0.2-mm accuracy.
Two homogeneous aquifer systems were simulated in this study, using glass beads of sizes 1,090 and 1,325 μm. The hydraulic conductivity of each of these was estimated using in situ measurement within the experimental flow tank. Various hydraulic gradients were successively imposed to the system, and the corresponding volumetric freshwater discharge was measured. The average hydraulic conductivity value was then subsequently derived using Darcy's law, as described in Oostrom, Hayworth, Dane, and Güven (1992). The average hydraulic conductivity value of the flow tank was estimated at 85 and 108 cm/min, for the beads 1,090 and 1,325 μm, respectively. The porosity of porous media was measured using volumetric method, and an average value of 0.38 was considered for all the beads.
The study included a total of 30 experiments, as shown in Table 1.
All the experiments were performed for the two bead sizes 1,090 and 1,325 μm. For each bead size, 10 experiments were completed for the

| Experimental procedure
A photograph of the whole experimental set-up used in this study is shown in Figure 2. The glass beads were packed into the central chamber under saturated conditions to avoid air bubbles development. The beads were packed in successive even-sized layers, and each layer was carefully compacted to ensure homogeneity of the packing. Two LED lights (Camtree 600) were placed behind the experimental set-up to illuminate the porous media, and a diffuser was fixed to the back of the tank to ensure homogeneity of the light throughout the aquifer.
A high-speed camera IDT MotionPro X-Series was used to capture images of the experiments every 30 s for all the experiments. Each time the camera was triggered, 10 images were recorded. The average of these images was then used in the analysis procedure.
To allow the determination the key SWI parameters, each experimental configuration was preceded by a calibration procedure allowing correlation of the saltwater concentration to the intensity of the light transmitting through the main chamber. This calibration is described in thorough details in Robinson et al. (2015). The purpose of the calibration was basically to correlate the saltwater concentration to the intensity of the light transmitting through the main chamber. The calibration procedure consisted in flushing the entire system with various saltwater solutions at known concentrations and recording the image of each saturated aquifer. For every concentration, the light intensity of every single pixel of the system was recorded. The known concentrations and corresponding light intensities were used to derive, through regression analysis, the coefficients a, b, and c of a power law curve expressed by C = aI b − c, adopted to relate the light intensity I and the concentration. The regression analysis was performed through a MATLAB code developed specifically for this purpose (Robinson et al., 2015).
The images recorded over the course of the experiments could thereafter be analysed through another MATLAB code that applied the derived regression coefficients determined in the calibration and analysed each single image to calculate the main intrusion parameters (i.e., toe length and width of TZ) and provided postprocessed images displaying the distribution of salt concentration throughout the system for clear visualization. The toe length was taken as the distance between the 50% concentration isoline and the saltwater boundary along the bottom boundary, and the width of the TZ was measured as the average of the vertical distance between the 25% and 75% saltwater isolines within the range 20% and 80% of the toe length.
Prior to the experiments, freshwater was injected at constant rate from a large tank located above the left side reservoir. The freshwater level was set high enough to allow the entire porous media to remain fully saturated with freshwater. Freshwater flux transited through the system from the inland boundary and exited at the coastal boundary without overflowing. The freshwater exiting the system was rapidly drained out through the overflow outlet of the saltwater reservoir.
The later was adjusted such that to maintain a constant saltwater head of 129.7 mm. Excess amount of saltwater was supplied from another large tank into the right reservoir to ensure the prompt flushing out of any freshwater floating at the surface. The density in the saltwater reservoir was continuously monitored using a hydrometer. FIGURE 2 Photograph of the experimental set-up: 1 = porous media chamber; 2 = freshwater reservoir; 3 = saltwater reservoir; 4 = ultrasonic sensors; 5 = high-speed camera; 6 = LED lights The saltwater intrusion experiments were initiated only after stabilization of the density measurements.
The experiments were initiated by lowering the overflow outlet of the freshwater reservoir such that to impose an initial constant freshwater head of 135.7 mm. This initial head boundary difference dh = 6 mm allowed the dense saline water to penetrate into the porous media until the system reached steady-state condition.
Steady-state condition was assumed when no substantial changes could be visually observed in the freshwater-saltwater interface position. Various hydraulic gradients were thereafter applied to the system to simulate changes in the hydrological conditions. This was done by varying the freshwater level such that to impose various head differences ranging from dh = 6 mm to dh = 3.6 mm, corresponding to hydraulic gradients of 0.0158 and 0.0095, respectively. These values were within the range of hydraulic gradients values commonly used in previous similar laboratory studies (Chang & Clement, 2012;Goswami & Clement, 2007) as well as within the range of values measured in some real coastal aquifers (Attanayake & Sholley, 2007;Ferguson & Gleeson, 2012).

| Description of the numerical model
The MODFLOW family variable-density flow code SEAWAT (Guo & Langevin, 2002) was used to perform the numerical modelling experiments. SEAWAT combines the modified MODFLOW (Harbaugh, Banta, Hill, & McDonald, 2000) and MT3DMS (Zheng & Wang, 1999) into a single program that solves the coupled groundwater flow and solute transport equations. SEAWAT has been widely used to solve various variable-density groundwater flow, such as the two box problems, the Henry problem, the Elder problem, and the HYDROCOIN problem (Guo & Langevin, 2002;Langevin, Shoemaker, & Guo, 2003). SEAWAT has also been successfully used for modelling SWI experiments involving freshwater head boundary variations (Chang & Clement, 2012;Goswami & Clement, 2007).
The purpose of performing numerical simulations was primarily to assess the consistency of the experimental data with the numerical predictions.
The  Abarca and Clement (2009). The spatial discretization satisfies the criterion of numerical stability; that is, grid Peclet number was less than or equal to four (Voss & Souza, 1987).
The porosity was set to 0.38. The molecular diffusion was neglected in all the numerical simulations (Riva, Guadagnini, & Dell'Oca, 2015).
The specific storage was set at 10 −6 cm −1 . The freshwater density was set to 1,000 kg/m 3 , and the saltwater density was set to 1025 kg/m 3 . Considering the standard density-concentration slope factor of 0.7, a concentration of 36.16 g/L was used for the seawater boundary.
3 | RESULTS AND DISCUSSION

| External SWI metrics
The experimental images of the steady-state saltwater wedge for the two cases are shown in Figure 3. In the first experiment, the head difference was varied as a function of the freshwater head boundary, whereas the head at the seaside boundary remained fixed at 129.7 mm. As stated earlier, the first steady-state saltwater wedge after the freshwater boundary head was dropped to 135.7 mm (dh = 6 mm) was considered as the initial condition of the experiment In the three following stress periods, the freshwater head was successively decreased to establish dh = 5.2, 4.4, and 3.6 mm in each of the three stress periods, respectively. In the last stress period, the freshwater was reset to the initial value (dh = 6 mm), which forced the retreat of the saltwater wedge towards the boundary. Figure 4 presents the comparison between the numerical simulation results and the experimental toe length data. In overall, the results show that the transient toe movement was reasonably well predicted by the numerical model. The maximum percentage difference was 10%, 1%, and 7% for dh = 5.2, 4.4, and 3.6 mm, respectively. Some mismatch could be observed in the retreat, where the toe movement predicted by the numerical was slightly faster than in the physical model, as observed in previous similar studies (Chang & Clement, 2012;Robinson et al., 2015;Robinson et al., 2016), with a maximum percentage difference recorded at the end of the retreat (nearly 21%).
The transient toe length data for the two cases are shown in The time to reach steady state corresponds to the time at which ΔTL becomes smaller than 1%. The curves of ΔTL of the intruding and receding wedge phases for the various water-level changes are shown in Figure 6 and the times required for the saltwater toe to reach steady state are presented in Table 2.   The data in Figure 6 show that the receding rate of the saltwater wedge was relatively faster than the intruding rate in both SL and FW scenarios, with an intruding wedge requiring up to twice longer time than the receding wedge to reach steady-state condition in case 1,325. Although previous studies demonstrated the asymmetry of timescale between intruding and receding processes following inland head variations (Chang & Clement, 2012;Robinson et al., 2015), the present results suggest that timescale asymmetry would also apply in scenarios involving sea-level fluctuations. Figure 6 also shows that the receding migration rate was relatively faster in case 1,325 than in case 1,090 in both SL and FW scenarios, whereas the intruding rates remain relatively comparable. This can readily be observed in the steepness of the slope at t = 0 min in case 1,325, as well as the shorter time taken to reach steady-state condition in both SL and FW fluctuation scenarios in Table 2. This means that for equivalent magnitude of freshwater-level rise/saltwater-level drop, faster seaward motion of the saltwater wedge would occur in higher hydraulic conductivity media, which is extension of the results of Robinson et al. (2016) to sea-level change scenarios. This trend was also confirmed in Figure 7 and Table 3, where the ΔTL results recorded for two different inland head change magnitudes but similar head difference are presented.

| Freshwater-saltwater TZ
The comparison between the transient experimental and numerical results is presented in Figure 8. Conversely, Figure 9 shows that the width of the TZ exhibited a noticeable widening during the saltwater retreat process, with a magnitude increasing with increasing freshwater-level increment. The peaks of the widening along with their time of occurrence are shown in Table 4. The results not only show that the widening process seems    TZ. The results also show that the TZ occurring during the retreat was also larger when the grain size was smaller (case 1,090). The data show that from the highest (case 1,325) to the lowest grain size (case 1,090), the peak of the widening increased by 33%, 21% and 43%, for ΔH = 0.8, 1.6, and 2.4 mm, respectively.
The response of the width of the TZ to freshwater-and saltwaterlevel fluctuations during the advancing wedge conditions was analysed, and the results are presented in Figure 10. This further highlights the key role of laboratory experiments in the investigation of density-driven flow in porous media and the evaluation of numerical models (Stoeckl et al., 2015).
It is hypothesized that this yellowish pulse observed in the experiments may have been the results of some freshwater "sliding" into the receding saltwater wedge along the bottom boundary, which was thereafter diluted and sucked upwards towards the interface due to density contrast, thereby exacerbating the widening. This possible sliding freshwater into the wedge may have been caused by the rather sudden rise of the freshwater head that induced in great change of the freshwater flow velocity in the system, thus abruptly pushing seaward the saltwater wedge and creating a perturbation allowing a tiny portion of freshwater to leak along the aquifer bottom. This experimental observation was not an artefact of any wall effects (i.e., 2D vs. 3D transport) given the relatively small width of the tank (10 mm). Despite field evidence of such mechanism have never been documented in previous studies, it is reasonable to hypothesize that such mechanism could conceivably occur in real coastal aquifer systems especially where the non-uniformity aquifer bedrocks may enhance preferential flow. Further investigations would nonetheless be required to examine the reproducibility of this mechanism within larger scale coastal aquifer models, which would enable easier examination of the TZ dynamics.

| SUMMARY AND CONCLUSIONS
This study provided a quantitative and qualitative analyses of the effect of boundary water-level variations on saltwater intrusion dynamics and the temporal response of freshwater-saltwater TZ.
Two homogeneous cases of different bead sizes (and hence different hydraulic conductivity) were analysed. The main findings of the study are as follows: • The expansion of the TZ did not occur during the intruding phase of saltwater wedge whether it resulted from freshwater-level drop or saltwater-level rise, regardless of the magnitude of the head change. This finding gives initial indications that in coastal aquifer systems, the widening of the TZ width may not be substantial • By contrast, a substantial widening of the TZ width was observed during the retreat process, with a peak magnitude and a peak time relatively similar in the scenario involving a rise of the freshwater level relative to the scenario involving a drop of the saltwater level, in cases of equivalent absolute head change magnitudes.
This finding suggests that in coastal aquifer systems, receding saltwater wedges induced by a rise of the water table would tend to exhibit TZ expansions of the same magnitude and the same timescale of occurrence as when following a drop of the sea level, provided that the head change magnitudes are comparable but opposite directions.
• The magnitude of the widening of the TZ observed during saltwater retreat was larger and extended over longer time frame in homogeneous aquifers with smaller hydraulic conductivity, irrespectively of whether the head change occurred at the freshwater or saltwater boundary.
• The experimental observations revealed that the widening mechanism was also enhanced by the presence of some freshwater sliding and into the wedge during the saltwater retreat process, which was thereafter sucked upward because of density difference effects and eventually travelled along the TZ towards the outlet, thereby promoting further widening. The observed phenomenon could not be reproduced by the numerical simulation, which underlines the importance of physical model experiments for the investigation of density-driven flow in porous media in improving general process understanding. It seems conceivable that such mechanism could possibly occur in real coastal aquifer systems where the common non-uniformity of aquifer bedrocks could further enhance preferential flow and thus facilitate the leakage of some freshwater below the wedge, albeit field evidence of such effect have never been reported in previous studies. Further investigations involving larger aquifer models would enable further examination of the TZ dynamics and explore the reproducibility of this observed phenomenon at larger scale.
• The transient analysis of the saltwater toe length revealed that in cases of equivalent saltwater-level variations but opposite directions, the intruding wedge required up to twice longer time to reach steady-state condition than the receding wedge. This observation shows that timescale asymmetry between saltwater intrusion and retreat processes would also occur in scenarios involving sea-level fluctuations. When comparing freshwater and saltwater change scenarios, the intruding and receding rates of the saltwater wedge as well as the time taken to reach steady state were fairly comparable.