The impact of aquifer stratification on saltwater intrusion characteristics. Comprehensive laboratory and numerical study

Laboratory experiments and numerical simulations were utilized in this study to assess the impact of aquifer stratification on saltwater intrusion. Three homogeneous and six layered aquifers were investigated. Image processing algorithms facilitated the precise calculation of saltwater wedge toe length, width of the mixing zone, and angle of intrusion. It was concluded that the length of intrusion in stratified aquifers is predominantly a function of permeability contrast, total aquifer transmissivity and the number of heterogeneous layers, being positively correlated to all three. When a lower permeability layer overlays or underlays more permeable zones its mixing zone widens, while it becomes thinner for the higher permeability strata. The change in the width of the mixing zone (WMZ) is positively correlated to permeability contrast, while it applies to all strata irrespectively of their relative vertical position in the aquifer. Variations in the applied hydraulic head causes the transient widening of WMZ. These peak WMZ values are larger during saltwater retreat and are negatively correlated to the layer's permeability and distance from the aquifer's bottom. Moreover, steeper angles of intrusion are observed in cases where low permeability layers overlay more permeable strata, and milder ones in the inverse aquifer setups. The presence of a low permeability upper layer results in the confinement of the saltwater wedge in the lower part of the stratified aquifer. This occurs until a critical hydraulic head difference is applied to the system. This hydraulic gradient value was found to be a function of layer width and permeability contrast alike.


| INTRODUCTION
Saltwater intrusion (SWI) is initiated because of the density difference between seawater and freshwater in coastal aquifers. Projected increases in water consumption within coastal communities alongside sea level rise will intensify the phenomenon. Therefore, the study of SWI constitutes a necessary step towards facilitating the protection of freshwater reserves. Fieldwork investigations of saltwater intrusion have highlighted the impact of heterogeneous structures on the extent of saltwater intrusion (Kazakis et al., 2016;Li et al., 2009;McInnis et al., 2013;Weinstein et al., 2007). Since real-life hydrological systems are relatively complex, their study alone is not ideal for deriving generalized conclusions about the impact of heterogeneity on saltwater intrusion characteristics. Approximations of the actual phenomenon, such as physical laboratory experiments, numerical modelling or analytical solutions have been successfully utilized instead.
The impact of aquifer heterogeneity on SWI has been investigated using predominantly analytical solutions (Ketabchi et al., 2014;Strack & Ausk, 2015) and numerical modelling (Held et al., 2005;Ward et al., 2008). Analytical solutions presented by Strack et al. (2016) indicated that a decrease in the permeability (k) of the upper part of an aquifer can significantly reduce the total length of intrusion. Expanding this study, Shi et al. (2018) conducted an analytical investigation of the impact of sea level rise on SWI in layered aquifers. This research concluded that in head-controlled systems toe length was a function of the relative hydraulic conductivities between layers. Rathore et al. (2018) introduced a new parameter, called transmissivity centroid elevation (TCE), to better evaluate the effect of stratification in SWI characteristics. TCE was defined as the summation of the product of the transmissivity and elevation (above the aquifer base) of each layer divided by the total aquifer transmissivity.
Utilizing an analytical solution, it was concluded that for the same values of aquifer transmissivity, higher TCE values corresponded to longer toe length in head-controlled systems. Lu et al. (2013) conducted a numerical evaluation of the impact of stratification in the values of steady-state width of the mixing zone (the dispersion zone where the salinity of water varies from zero to 100%). Based on a setup of three-layered aquifers, it demonstrated that when a low permeability layer overlies a relatively high-k layer, the mixing zone in the low-k layer is widened. Mixing zone gets thinner in cases where highly permeable layers overlie less permeable ones.
To the best of our knowledge, there has never been, a laboratorybased study quantifying the impact of stratification on all three fundamental saltwater intrusion characteristics, namely the toe length (TL), width of the mixing zone (WMZ) and the angle of intrusion (AOI).
Although a few experimental studies were presented in the literature, those were based on steady state conditions (Lu et al., 2013) or were focusing on the effect of cut-off wall barriers (Abdoulhalik & Ahmed, 2017a). Unlike previous experimental studies that relied on simple visual observation of the sandbox data, the present investigation utilized the automated image processing techniques developed in Queen's University Belfast (Etsias, Hamill, Benner, Aguila, McDonnell, Ahmed, & Flynn, 2020a;Etsias, Hamill, Benner, Aguila, McDonnell, & Flynn, 2020b;Robinson et al., 2015). This has enabled the accurate reproduction of laboratory saltwater concentration fields, as well as the quantification of specific saltwater intrusion characteristics (TL, WMZ and AOI) with high accuracy. All the laboratory results were successfully recreated using numerical simulations. The trends identified based on the analysis of the sandbox data, were further expanded using numerical simulations, through six sensitivity analysis scenarios.
The main novel aspects of this study can be summarized as follows: • With a total of nine investigated aquifers, this is, to date, the most comprehensive laboratory study on the impact of stratification on saltwater dynamics. The direct comparison of saltwater hydrodynamics in different heterogeneous setups, enabled the identification of a series of new conclusions that have never been presented before.
• This is the only sandbox study incorporating high precision quantification of the experimental WMZ into each specific permeability layer of stratified aquifers. Supplementary to that, it assessed, for F I G U R E 1 3D representation of the utilized experimental setup the first time ever in a heterogeneous aquifer setup, the variation of WMZ under transient conditions.
• The impact of aquifer stratification on the AOI, was investigated for the first time in this study, under both steady state and transient conditions.
• The numerical sensitivity analysis identified the impact of four distinct aquifer characteristics: permeability contrast between strata; total aquifer transmissivity; total number and width of individual permeability layers. The effect of those four variables on saltwater intrusion dynamics was studied in isolation from one F I G U R E 2 (i) Photographs and the equivalent (ii) heterogeneous structure fields (generated by the classification artificial neural network) of the six investigated stratified aquifers another using numerical modelling. This is an element that was not taken into account in previous numerical investigations and has enabled the further expansion of the currently accepted concepts of SWI.
The conclusions of this study, derived for idealized twodimensional laboratory experiments, should be a valuable contribution towards better understanding of SWI dynamics in real life layered aquifer systems.

| LABORATORY SETUP AND TEST CASES
The laboratory setup utilized to derive the experimental data of the current investigation was described in detail by Robinson et al. (2015).
The sandbox apparatus consisted of a 0.38 m × 0.15 m × 0.01 m viewing chamber flanked by two cylindrical reservoirs ( Figure 1). The left tank was filled with freshwater, while dyed saltwater filled the right one. Head level in the side reservoirs was monitored by two ultrasonic sensors, while adjustable outflows enabled the application of different hydraulic gradients. Glass beads of three different diameters, 780, 1090, and 1325 μm, were siphoned into the viewing chamber to minimize the amount of air trapped inside the porous medium and ensure saturated conditions. Acrylic mesh screens stabilized the beads into the central chamber. Experiments were conducted in a dark room where illumination was provided by two LED light panels positioned on the back of the laboratory rig. High resolution multicolour images were acquired every 5 min with a Nikon D850 DSLR Camera.
Nine experimental aquifers were investigated, including three homogeneous aquifers, one for every bead size used, and six layered heterogeneous ones ( Figure 2). The setup of the stratified aquifers enabled the study of various degrees of heterogeneity.
The investigated layered aquifers can be divided into three categories: • Two-layered aquifers: in the first a low permeability layer overlaid a higher permeability one (780-1325 μm), while in the second the reverse setup applied where 1325 μm high permeability stratum was laid on top of the 780 μm low-k layer (Figure 2 (a),(b)).
In all tested aquifers the various strata were formed by glass beads of similar size. The layers created by the 780 μm glass beads were characterized by low (L) permeability, the 1090 μm layers by medium permeability (M) and the 1325 μm ones by high (H) permeability. In the interests of clarity, the six laboratory heterogeneous aquifers will be henceforth referred to according to the quali-  Initially the aquifer contained only freshwater. Three hydraulic head level differences (dH) were applied between the fresh and saltwater tanks. A head difference of 6 mm was applied to the system triggering the initial stage of saltwater intrusion. Subsequently, the head difference decreased (4 mm) resulting in a new phase of intrusion. Finally, saltwater retreat was initiated by an increase of the hydraulic head difference (5 mm). The hydraulic gradients investigated here are within the range of hydraulic gradients observed in some real world problems (Attanayake & Sholley, 2007;Ferguson & Gleeson, 2012). During the experiment, water in the right cylinder was monitored using a YSI Professional Plus Instrument (Pro Plus) water quality meter; this ensured that salinity did not reduce due to the outflow of freshwater through the porous medium. Total data acquisition time per intrusion or retreat phase varied from 50 to 90 min (Table 1), depending on the aquifers' heterogeneous structure.
Experimental images were post-processed using the machine learning oriented method introduced by Etsias, Hamill, Benner, Aguila, McDonnell, and Flynn (2020a). A multi-layered classification artificial neural network (ANN) was utilized to precisely recognize the heterogeneous structure of the six stratified aquifers (Figure 2). Subsequently, saltwater concentration fields were generated for all the steady-state saltwater intrusion images (dH = 6, 4 and 5 mm) by a 10-neuron regression ANN with one hidden layer ( Figure 3). The values of specific saltwater intrusion variables such as toe length of intruding wedge, width of mixing zone and angle of intrusion were calculated for the derived concentration fields by applying the data processing algorithms presented by Robinson et al. (2015). TL quantifies the extent of saltwater intrusion, which was determined as the distance of the 50% saltwater concentration isoline from the right (saltwater) edge of the laboratory aquifer. WMZ equalled to the average vertical distance between the 25% and 75% salinity isolines, while the AOI is the angle that the saline wedge formed with the bottom of the aquifer, or with the interface of different permeability layers.  (Voss & Souza, 1987). Hydrostatic pressure was specified at the right (saltwater) and left (freshwater) boundaries of the model to simulate the head level differences of 6, 4 and 5 mm. An intrinsic flow test on the experimental domain allowed the calculation of the permeability of the porous media using Darcy's law. Saturated-unsaturated groundwater flow models were performed since similar models have been successfully utilized before to simulate saltwater intrusion occurring in a sandbox setup (Houben et al., 2018).
In this study, the van Genuchten equation (van Genuchten, 1980) was used to simulate the unsaturated flow. The van Genuchten parameters of the numerical model were derived from experimental data for glass beads of comparable size (Benson et al., 2014;Sweijen et al., 2017). The values of applied dispersivities conformed to the range determined by Abarca and Clement (2009). Model input parameters are summarized in Table 2. The simulation time for each recreated laboratory case was identical to the equivalent experimental duration, while the time-step was equal to 1 s.
The study was expanded via sensitivity analysis that evolved around six simulation scenarios. The total simulation time for each intrusion and retreat phase on these scenarios was 60 min. During this period, steady state was achieved in all numerical cases. The layer permeability values used here were kept similar to those measured in the laboratory setup. In none of the simulations did the utilized permeability values vary more than an order of magnitude from the experimental ones. To better visualize the correlation between the SWI characteristics and aquifer stratification variables, nonlinear fitting was applied between the numerical results in all scenarios.
Since formulating precise generalized equations was beyond the scope of this study, power law was the basis for the conducted interpolations, and the data fit was deemed satisfactory in all cases. For the sake of clarity, these scenarios are presented alongside the experimental and numerical results in Section 4. The heterogeneous structure of the layered aquifers investigated during the six sensitivity analysis scenarios is presented in detail in the Appendix S1.

| Experimental results
The experimental steady-state saltwater concentration fields of the nine investigated aquifers are presented in Figure 3. The laboratory data were successfully reproduced with numerical simulations, during both the two intrusion and the final saltwater retreat phase (Figures 4 and 5 and A1 of the appendix). In contrast to the homogeneous aquifer settings, experimental saltwater wedges in the layered cases demonstrated significant refraction at the freshwater-saltwater interface.
The general shape of the simulated saltwater wedges as well as the specific saltwater characteristics that is, TL, WMZ and AOI exhibited a high quality fit with the experimental data, nevertheless some variation could still be identified. The TL values generated by the utilized numerical model ( Figure 5) were comparable to the experimental TLs for the initial intrusion (dH = 6 mm) and the final retreat phase (dH = 5 mm), presenting an average relative difference of 5.3% and 4.8% respectively. During the second intrusion phase (dH = 4 mm) the numerical model underpredicted the extent of saline intrusion by an average of 9.7%, while the difference between experimental and numerical results was the greatest for the H-M-L aquifer (5.6 cm).
Similarly, the numerical saltwater concentration fields in the H-L aquifer ( Figure 4(b)) included a far thinner mixing zone inside the low permeability stratum than its equivalent experimental field.
In sandbox setups of relatively small size, experimental errors can be introduced by flaws in water level adjustment and sandbox levelling, as well as by uneven flow through the side mesh screens or variations in the density of the saline solution. Moreover, minor heterogeneities due to small variations in the diameter of glass beads and packing may exist in the experimental settings while not accounted for in the numerical model, which assumes perfect homogeneity for each bead layer. Automated image analysis algorithms were employed in this study to determine the boundaries between the layers. Nevertheless, since the size of utilized elements was similar to that of a single glass bead, there was some divergence between the generated numerical and experimental permeability fields. Finally, numerical TL underprediction, of comparable scale, has been reported in previous The simulated velocity vector fields for the six heterogeneous laboratory aquifers are presented in Figure 6. In addition to validating The acquired experimental data allowed novel concepts about the impact of stratification on SWI to be formulated and subsequently expanded through rigorous sensitivity analysis. The effect of stratification on each one of the main SWI characteristics is now presented.

| Toe length
For all experimental images, the toe length was inversely correlated to the applied hydraulic head difference, which is in agreement with previous investigations (Houben et al., 2018). The lengths of the experi- to quantify the impact of stratification on SWI. TCE values alongside average aquifer permeability were calculated for the six heterogeneous experimental aquifers (Table 3) using the following equations: where n is the number of layers, k is permeability, b is layer width and y is layer elevation above the aquifer's base.
The laboratory data indicated that aquifers with higher TCE and average permeabilities manifested longer saline wedges (Table 3).
Even though TCE was formulated using idealized steady state analytical solutions, which assumed a sharp interface between freshwa-  Table 4. It is evident that in all the heterogeneous cases, WMZ of the L layers was wider than that of the equivalent L homogeneous aquifer, while for the H bead size WMZ decreased. These data expand the conclusion of Lu et al. (2013) that identified this variation of WMZ strictly on the upper-   Finally, a clear correlation was established between the layer's relative position, and the time in which the peak WMZ was observed. In all six stratified aquifers, widening occurred first in the lower strata, while significant time lagging, was observed for the layers closer to the aquifer's free surface. This time lagging, between the change of the hydraulic gradient and the corresponding peak WMZ, varied significantly from aquifer to aquifer, ranging from 15 up to 55 minutes.
The observed time lagging in each hydrological system was in agreement with the total experimental time that was deemed necessary for the stabilization of the saline wedge (Table 1). In accordance with that, the study, is necessary. The current study of AOI constitutes a contribution towards this direction.
In general, the laboratory data demonstrated that saltwaterfreshwater interface slope is mild in high permeability zones and steep in low permeability zones, this is in agreement with the numerical and laboratory investigations of Abarca (2006) and Robinson et al. (2016).
The AOI values for the quasi steady-state conditions under dH = 4 mm, are presented in Table 6. The steepest angle of intrusion,   A1f and A1h of the appendix). It was observed that given a specific k contrast between two layers, saltwater would not enter into the upper layer until a certain value of hydraulic head difference was applied to the system. This value, identified for the first time in this investigation, will be henceforth called critical head difference (dH crit ).
When a greater head difference (dH = 6 mm) is applied to the hydrological system, more freshwater flows into and subsequently out of the aquifer. Since in these L-H setups permeability is limited in the upper part of the aquifer, the outflow area widens incorporating parts of the more permeable lower zone. This confines the saline wedge inside the lower stratum. Only in cases when the inflow gets reduced to a sufficient level, can outflow occur strictly through the upper layer.
At this point, saltwater intrudes in the upper portion of the aquifer.
These results could be important for multiple real-life applications. Intrusion of saltwater on the upper-most parts of freshwater aquifers negatively impacts soil quality thus affecting human activities such as farming (Alam et al., 2017). Moreover, one of the main variables influencing pumping induced saline intrusion is the absolute distance between the pumping wells and the saltwater-freshwater interface (Abdelgawad et al., 2018;Abdoulhalik & Ahmed, 2018b).
Thus, being able to safely predict the pre-pumping distance of the saline front from the aquifer's free surface is a prerequisite for the successful management of actual coastal aquifers. The correlation between dH crit and layer width and permeability was further quantified using numerical simulations (Scenario 6).

| Sensitivity analysis
4.2.1 | Sensitivity analysis Scenario 1 (effect of permeability contrast on the TL) Six two-layered numerical models were created with varying k contrast values between their layers (Table A1 of the appendix). In half of these settings, a high permeability layer overlaid a lower permeability one (H-L cases), while in the remaining setups the upper layer permeability was lower (L-H cases). All simulated heterogeneous aquifers had the same average permeability, 2.39 × 10 −9 m 2 , which corresponded to the permeability of 1325 μm beads homogeneous aquifer. Four distinct head difference values (dH = 5, 6, 7 and 8 mm) were applied on the models. The calculated toe lengths are presented in Figure 11.
It was established that for equal total aquifer transmissivity, higher k contrast (k upper /k lower ) corresponds to more extended saltwater intrusion while less permeable zones in the upper part of the aquifers drive more freshwater towards lower layers ( Figure 6). The presence of more freshwater near the aquifers' bottom contributes to the seaward repulsion of the saline wedge. The smaller the value of k contrast , as is the case in the low -high settings, the more intense is the containment of the saline wedge and vice versa. The toe length was negatively correlated to the applied head difference, while there was a linear relationship between toe length and the natural logarithm of k contrast . The slopes in the linear equations of Figure 11(b) express the F I G U R E 1 1 Impact of layer permeability contrast on the steady state toe length of saltwater wedge on simulated two-layered heterogeneous aquifers with the same average aquifer permeability (2.39 × 10 −9 m 2 ). Toe length plotted against (a) head difference and (b) logarithm of k contrast impact of k contrast to TL, while the second term represents the influence of the applied head difference. It is obvious that for smaller head differences the values of the slope tend to plateau (around 0.018 for the investigated aquifers) while the constant component gradually increases.

| Sensitivity analysis Scenario 2 (effect of total aquifer transmissivity on the TL)
Five two-layered numerical models were developed with varying average aquifer permeability but maintaining the same permeability contrast between their layers (k upper /k lower = 2) ( Table A2 of  The evaluation of experimental data, alongside a detailed sensitivity analysis, indicated that aquifer stratification significantly affects the length of intruding wedge; this is in agreement with previous studies (Abdoulhalik & Ahmed, 2017a;Lu et al., 2013). The study expanded the current knowledge, proving that the extend of saline intrusion in stratified systems is a multivariable problem, more complex than what previously hypothesized. In particular, TL in stratified systems is positively correlated to three different aquifer characteristics: the permeability contrast between layers, the total aquifer permeability and the number of different strata. Taking into account the effect of all three aforementioned stratification characteristics is crucial in the successful evaluation of the physical extent of the intruding saltwater wedge length.

| Sensitivity analysis Scenario 4 (effect of permeability contrast on the steady-state WMZ)
This scenario was based on two-layered numerical models (Table A3 of the appendix). In the first model (Figure 14 4.2.6 | Sensitivity analysis Scenario 6 (effect of layer width and permeability on the critical head) Three setups of two-layered aquifers were simulated. In all cases a low permeability layer overlaid a higher permeability one. Upper layer permeability was kept constant equal to 7.89 × 10 −10 m 2 while permeability of the lower layer was 1.5-10 times larger (1.19-7.89 × 10 −9 ). In each setup, the width of the upper layer varied such that it occupied 3 /4, ½ and 1 /4 of the total aquifer width. The layer permeability values of these numerical setups are presented in Table A5 of the appendix. It was found that the critical head difference was negatively correlated to k contrast in a similar way for all three cases ( Figure 16). The width of the upper layer induced an important impact on both first and third coefficient of the power law equation. This is the first study that defined this critical head difference and correlated its values with aquifer heterogeneity characteristics, such as the width of the different permeability strata and their permeability contrast.

Toe length
Three heterogeneity characteristics were identified as the ones with dominant effect on the length of intrusion, including the permeability contrast between layers, total aquifer transmissivity and the number of layers. The following conclusions were drawn: • Heterogeneous aquifers with more permeable upper layers (H-L cases) exhibited the longest saltwater intrusion amongst layered aquifers with similar average permeability. A positive linear relationship correlated the logarithm of k contrast with the generated TL.
• In stratified aquifers with the same permeability contrast between layers, higher total aquifer transmissivity induced longer saltwater intrusion. The larger was the applied head difference, the smaller was the impact of transmissivity.
• In layered aquifers with the same k contrast and total aquifer transmissivity, the toe length was positively correlated to the number of permeability layers via a power law function. The upper limit of this F I G U R E 1 6 Impact of layer permeability contrast on the critical head difference for three distinct upper layer width setups, where the upper layer occupies 1 /4 (red), ½ (blue) or 3 /4 (black) of the total aquifer depth function was defined by the TL in systems were aquifer permeability changed linearly with depth.
These findings expand the established concept that TL in headcontrolled stratified systems is a mere function of permeability contrast between strata (Shi et al., 2018). Instead, its relationship with aquifer stratification characteristics was proven to be much more complex than what was previously hypothesized.

Width of the mixing zone
• Expanding already established concepts about steady state WMZ, it was proven that when a lower permeability layer was either overlying or underlying higher permeability ones, its mixing zone was widened. On the other hand, WMZ was reduced in the case of the more permeable layer. The change in WMZ was positively correlated to k contrast . Even though the identified WMZ variation was substantially greater for the overlying layer, the investigation proved that, this variation applies for all strata irrespectively of their relative vertical position in the aquifer, further expanding on the current concepts.
• Every change on the applied hydraulic head difference, was followed by an explicit transient widening of the mixing zone. This peak WMZ was considerably larger during the saltwater retreat phase of the experiments. Peak WMZ values were negatively correlated to each stratum's permeability, and to its vertical distance from the aquifers' impermeable bottom. Finally, the widening of the mixing zone was not instantaneous, instead a distinct time span was identified between the dH re-adjustment and the peak WMZ values. This time lagging was considerably longer for the uppermost permeability layers and shorter for the lower strata.

Angle of intrusion
• Experimental data revealed that the contrast between various permeability layers significantly altered the steady state values of AOI, especially in the upper part of the aquifers. In H-L setups, the AOI of the upper layer became smaller, while it increased in the inverse L-H cases. The relative change in AOI values was proportional to the permeability contrast between layers. Its effect being more significant in the L-H cases. This was the first time that the impact of stratification on AOI was quantified.
• At the point of initial intrusion of the saline wedge into an individual permeability layer, the transient saltwater-freshwater interface was up to 60% steeper than the interface's orientation after steady state was achieved. These steep AOI values, were negatively correlated to the layer's permeability, while they were observed only in the initial experimental phase and did not characterize any subsequent variations of the applied hydraulic gradient.

Critical head difference
• In cases were a low permeability layer overlaid more permeable aquifer structures, the saltwater wedge did not enter the top, less permeable layer. On the contrary it stabilized in its proximity, until a specific critical value of hydraulic head difference was applied to the system. This phenomenon, observed in three of the experimental stratified aquifers, was further investigated via numerical simulations. It was derived that the relationship between permeability contrast and critical head difference is an exponential formula, while the critical dH was negatively correlated to the width of the low permeability layer. Identifying this critical head difference and quantifying its relationship with stratification characteristics is the final novel contribution of the current study.
The current investigation outlined the impact of aquifer stratification on saltwater intrusion characteristics. Even though derived from an idealized two-dimensional laboratory setup, the novel conclusions of this study provide valuable insight on how realistic aquifer stratification could impact the main SWI characteristics, thus contributing towards the successful mitigation of the phenomenon. Although actual relationships between aquifer characteristics and the values of TL, WMZ and AOI were obtained via non-linear interpolation between the acquired numerical results, formulating generalized equations applicable in a wide variety of aquifer setups was beyond the scope of this work. Nevertheless, the presented findings could act as a solid foundation for future analytical studies in this direction. Finally, the detailed experimental data presented in this paper can be utilized in benchmarking models in upcoming numerical investigations of SWI in stratified coastal aquifers.