Flow and Transport in Coastal Aquifer‐Aquitard Systems: Experimental and Numerical Analysis

Coastal aquifers are commonly layered, and thus, a clear understanding of groundwater flow and salt transport in layered coastal aquifers is important for managing fresh groundwater. However, the influence of leakage between adjacent aquifers on flow and transport processes remains largely unknown where the influence of tides is considered. This study used laboratory experiments and numerical simulation to examine the processes of flow and transport within a tidal aquifer‐aquitard system (i.e., an unconfined aquifer underlain by a semi‐confined aquifer, with an intervening thin aquitard). The laboratory‐scale observations of the current study are the first observations of offshore fresh groundwater within a semi‐confined coastal aquifer. The numerical and laboratory results are in close agreement, revealing that upward leakage from the semi‐confined aquifer into the saltwater wedge of the overlying unconfined aquifer caused buoyant instabilities to form. The development of freshwater fingers created complex saltwater‐freshwater mixing, leading to mixed saltwater influx‐efflux patterns across the sloping aquifer‐ocean interface. Compared with non‐tidal conditions, tidal forces reduced the net upward leakage from the semi‐confined aquifer to the overlying unconfined aquifer. This increased the horizontal flow toward the sea, which in turn reduced the extent of the saltwater wedge in the semi‐confined aquifer. The higher rates of both fresh and saline submarine groundwater discharge (SGD), caused by tides, led to lower groundwater ages in the semi‐confined aquifer. These findings have important implications for unveiling the complex characteristics of seawater intrusion, SGD and geochemical hotspots within layered coastal aquifers.


Introduction
Coastal aquifers connect terrestrial and oceanic systems, and contain complex and dynamic hydrogeological processes driven by water density variations, tidal effects and time-varying hydrological forces from land-based stresses, such as recharge and pumping (Michael et al., 2017;Werner et al., 2013).Mixing and reaction processes occur between freshwater, seawater and other chemical constituents (e.g., nutrients, carbon, metals, etc.) within coastal aquifers, and these influence groundwater-based chemical loads to the sea.Additionally, tides and aquifer heterogeneities create complex pathways of coastal groundwater flow, with profound effects on groundwater residence times, groundwater-borne chemical loads, and the location of groundwater outflows (Geng et al., 2020;Heiss et al., 2017;Kim et al., 2020;Li et al., 1999).A comprehensive understanding of these (and other) flow and transport processes within coastal aquifers is a precursor to the sustainable management of fresh groundwater resources and the ecosystems of coastal zones.
The outflow of groundwater to the sea, known as submarine groundwater discharge (SGD), includes recirculating seawater and fresh groundwater that derives from recharge to inland catchments (Burnett et al., 2003;Michael et al., 2005;Moore et al., 2008).Tides may significantly enhance this seawater circulation, such that SGD (i.e., the sum of freshwater and saline groundwater discharge) may be substantially larger where tides are taken into account (Robinson, Li, & Barry, 2007).Tides also influence the mixing dynamics, including between freshwater and seawater, but also in regards to the fate of coastal zone contaminants from agricultural, industrial and urban sources (Heiss & Michael, 2014;Robinson, Gibbes, et al., 2007).
The impact of tides on seawater intrusion (i.e., the seawater extent in coastal aquifers) has also been widely investigated.Ataie-Ashtiani et al. (1999) were the first to show changes in the extent of seawater in unconfined coastal aquifers due to tides.Their numerical modeling results indicated that tides caused diluted seawater to intrude further inland compared to non-tidal conditions.However, subsequent studies of the effects of tides on the seawater extent within coastal aquifers have discovered the opposite effect, whereby tides are more likely to reduce the distance inland that saltwater occurs within the unconfined aquifer.This is attributable to the occurrence of an additional, tide-induced plume of saltwater, termed the "upper saline plume" (USP) by Robinson et al. (2006).The formation of the USP causes fresh groundwater to discharge to the shoreline at a lower elevation, exiting at around the low-tide mark (rather than at mean sea level, as occurs without tides).This is the primary mechanism for the saltwater wedge to shrink to a reduced landward extent (relative to the non-tidal case) in unconfined coastal aquifers (Fang et al., 2021;Kuan et al., 2012).Studies that have investigated the effects of tides on seawater intrusion in confined aquifers include the investigation of Inouchi et al. (1990).They used numerical modeling to reveal that the saltwater extent fluctuated with the tidal cycle, intruding furthest inland when the tidal sea level reached mean sea level during the falling tidal period.Chen and Hsu (2004) also used numerical modeling to study the effects of tides on coastal confined aquifers, finding that tides induced only minor fluctuations in the toe position, and that the time-averaged location of the saltwater wedge toe was largely insensitive to tidal effects.
Real-world observations of the USP have been obtained for unconfined coastal aquifers from around the world, including Moreton Island, Australia (Robinson et al., 2006), Waquoit Bay, USA (Abarca et al., 2013), Cape Henlopen, USA (Geng et al., 2017;Heiss & Michael, 2014) and Shilaoren Beach, China (Zhang, Wu, et al., 2021).The USP contains higher velocities than those of the density-driven saltwater circulation cell (Xin et al., 2010), that is, the saltwater wedge, which is contained in the lower part of the aquifer due to its higher density relative to freshwater.This results in higher levels of dissolved oxygen in the USP compared to those found in the saltwater wedge because the travel times are shorter for the tide-driven circulation (Robinson, Gibbes, et al., 2007).In general, the higher rates of flow within the USP create a more active area for geochemical reactions than the saltwater wedge (Robinson et al., 2009).
Homogeneous coastal aquifers are often assumed in conceptual models when exploring fundamental mechanisms of flow and transport in coastal aquifers subject to tidal fluctuations (Li et al., 1999;Robinson et al., 2009).In coastal regions containing layered aquifer systems, deeper aquifers may be semi-confined, creating the opportunity for inter-aquifer fluxes to occur (Kim et al., 2006;Love et al., 2000).Tides may cause fluctuating interlayer fluxes that modify the classical salinity patterns that are illustrated in single-aquifer investigations.While the importance of tidal forcing acting on coastal unconfined aquifers has been widely illustrated, few studies have explored the tidal influence on flow and salt transport in unconfined aquifers where there is upward leakage through the base.
Investigations of tidal effects in semi-confined aquifers include analytical solutions (Jiao & Tang, 1999;Li & Jiao, 2001), numerical modeling (Evans et al., 2020) and field investigations (Evans et al., 2020;Ratner-Narovlansky et al., 2020).Early studies of tidal, semi-confined aquifers focused primarily on hydraulic processes, such as the propagation of tides, and the occurrence of seawater in the aquifer was neglected.Jiao and Tang (1999) found that downward leakage from overlying unconfined aquifers considerably reduced the amplitude of tide-induced groundwater-head fluctuations within semi-confined aquifers, leading to smaller distances of tidal propagation.Their work was extended by Li and Jiao (2001) to account for the offshore extension of the aquifer that is common for deep, semi-confined coastal aquifers.They developed analytical solutions for the propagation of tides emanating from offshore, and revealed that tide-induced fluctuations in the onshore part of the semi-confined aquifer significantly decrease with overlying aquitards that extend further offshore.2020) conducted a regional-scale modeling study of flow and salt transport in a multi-layered coastal aquifer at North Inlet, South Carolina, USA.Their numerical results found that tidal fluctuations have a negligible influence on the flow and salinity distribution in offshore semi-confined aquifers far from the shoreline (i.e., 10-15 km), at least for the relatively low upward leakage through the overlying aquitard (i.e., leakage of approximately 1 mm yr 1 ) that they considered.The low rate of inter-aquifer leakage creates a rather weak hydraulic connection between two overlying aquifers.This likely restricted the downward propagation of tidal effects into the semi-confined aquifer.Where aquitards have higher hydraulic conductivity (K s ), as observed in other studies of layered-aquifer systems (e.g., Atteia et al., 2005;Knight et al., 2018;Love et al., 2000), larger leakage rates than those considered by Evans et al. (2020) are expected.How the salinity dynamics, including the seawater circulation and the salinity distribution, in semi-confined, coastal aquifers respond to aquitards of different K s , and with alternative amplitudes for tides, is not clear.To date, there have been no physical observations, such as those that can be obtained using laboratory experiments, to examine the tidal effects on the characteristics of seawater intrusion in semi-confined aquifers where there is upward leakage through the less permeable roof.Additionally, the temporal and spatial variability in inter-aquifer leakage induced by tidal fluctuations has not been previously assessed.These factors are likely important in controlling flow and transport within semi-confined coastal aquifer systems.

Evans et al. (
Figure 1 illustrates the conceptual model adopted in this study, comprising a two-dimensional aquifer-aquitard system that consists of an unconfined aquifer, an aquitard and a semi-confined aquifer.Here, both the aquitard and the semi-confined aquifer extend beneath the sea, and the unconfined aquifer is bounded by the sloping beach surface.The effects of inter-aquifer leakage and tides on the flow and transport within this system are evaluated in the current study using laboratory experiments and numerical simulation.The principle aims are to: (a) examine the influence of any tide-induced, inter-aquifer leakage on the saltwater distribution within the layered-aquifer system; (b) investigate the spatiotemporal distribution of submarine groundwater discharge to the seaward boundary, including for both non-tidal and tidal conditions; (c) assess the effect of tidal fluctuations on groundwater transit times and ages within both the unconfined and semi-confined aquifers.The relationship between salinity distributions and tidal dynamics is considered for layered aquifers with stronger hydraulic connections than those considered by Evans et al. (2020), consistent with the conditions found within many sites around the globe.

Laboratory Experiments
Laboratory experiments were conducted using a 4 m (length) × 0.02 m (width) × 0.8 m (height) sand flume, as shown in Figure 2a.The sand flume has three parts: the left compartment is a saltwater reservoir representing the seaward boundary, the middle compartment is packed with quartz sand, and the right compartment is a freshwater reservoir that represents the landward boundary.Two 0.01 m-thick stainless steel screens prevent the sand in the middle compartment from entering the two reservoirs.The narrow width of the sand tank creates a 2-D crosssection (perpendicular to the shoreline).To minimize the sidewall flow, the internal sidewalls of the tank were roughened, as done in previous studies (e.g., Lu et al., 2013).
The saltwater used in laboratory experiments was prepared by dissolving 36 kg of sodium chloride (NaCl) into 1 m 3 of deionized water at 20°C.The saltwater and freshwater (i.e., deionized water) densities were determined by a densimeter to be 1,024 ± 1 kg m 3 and 998 ± 1 kg m 3 , respectively.Red food dye (i.e., Allura Red AC, a red azo dye) was added to the saltwater for visualization.Three types of sand were used in the experiments: wellsorted coarse sand (i.e., unconfined and semi-confined aquifers), well-sorted medium sand (i.e., high-leakage aquitard), and well-sorted fine sand (i.e., low-leakage aquitard).Grain-size distribution curves for the three sand types produced d 50 = 1.76 mm and d 60 /d 10 = 1.88 for the coarse sand, d 50 = 0.35 mm and d 60 /d 10 = 1.74 for the medium sand, and d 50 = 0.09 mm and d 60 /d 10 = 1.65 for the fine sand.Here, d 10 , d 50 and d 60 are the grain sizes that 10%, 50% and 60% of samples are finer by weight, respectively.The saturated hydraulic conductivity (K s ) of the coarse, medium and fine sands was 1.94 × 10 2 m s 1 , 5.73 × 10 4 m s 1 and 3.89 × 10 5 m s 1 , respectively, obtained using in-situ, sand tank testing and application of Darcy's Law based on steady freshwater flows.The porosity (θ) of the coarse, medium and fine sands was obtained from the oven-drying method, which produced values of 0.42 ± 0.01, 0.45 ± 0.01, and 0.46 ± 0.01, respectively.
The landward reservoir was supplied with freshwater at a fixed rate of 64 mL min 1 through a peristaltic pump during all experiments.This inland constant flux was determined from a numerical simulation of the homogeneous coarse sand aquifer, with the inland boundary set to a constant head of 0.565 m, and in the absence of tides.The use of a constant-flux boundary meant that changes to the hydraulic conditions at the ocean boundary, such as the occurrence of tide-induced watertable overheight (i.e., the super-elevation of head conditions at the shoreline arising from tidal effects; see Nielsen, 1990), did not change the time-averaged freshwater component of the total SGD.
The water level at the seaward boundary was controlled by a variable-height overflow column, which was connected to a tidal generator in experiments requiring tides.In all cases, mean sea level was 0.54 m (from the internal base of the sand tank).Sinusoidal tides were created by the variable-height overflow with a period of 60 s and an amplitude of 0.041 m.Freshwater discharging from the aquifers into the saltwater reservoir diluted the saltwater salinity during the experiments.To maintain a relatively constant saltwater salinity within the reservoir, another peristaltic pump was used to skim the low-salinity water from the top of the saltwater reservoir, which was then discarded through a drain.
The sand aquifer was divided into three parts, including an offshore part (i.e., x = 0-0.6 m), a sloping beach (i.e., x = 0.6-1.6 m) and an onshore part (i.e., x = 1.6-4.0m), where x = 0 represents the offshore limit of the conceptual model (Figure 1).The tank was filled with sand to a depth of 0.34 m for the offshore part and 0.64 m for the onshore part, leading to a shoreline with a slope of 0.3.The top and bottom elevations (taken from the internal base of the sand tank) of the aquitard were 0.34 and 0.3 m, respectively.Three types of layered-aquifer settings were considered in the experiments: (a) a homogeneous aquifer (Case H, filled with coarse sand); (b) an aquiferaquitard system with the aquitard comprised of medium sand and a coarse sand aquifer (Case L1); (c) an aquiferaquitard system with fine sand as the aquitard and a coarse sand aquifer (Case L2).The naming convention for cases defines whether experiments were non-tidal (e.g., Cases H-S, L1-S, and L2-S, with "S" meaning static sea level) or tidal (e.g., Cases H-T, L1-T, and L2-T, with "T" indicating tides).Field-scale scenarios are Case F-S for the non-tidal condition and Case F-T for the tidal condition (see Section 2.2.2).Each laboratory experiment was continued until the measured variables (e.g., the toe position of the saltwater wedge) indicated that the system had reached steady-state (i.e., non-tidal conditions prior to the tide being introduced) or quasi-steady conditions (after tides are introduced).The experiments were recorded using a high-resolution digital camera (Canon IXUS 175) at a time interval of 5 min.The indicator for the steady and quasi-steady states in the laboratory experiments is the toe position of the saltwater wedge remaining unchanged within 5 hr.

Model Setup
SUTRA-MS (Hughes & Sanford, 2005), a variable-saturation, variable-density groundwater flow and multiplespecies solute transport model, was used to simulate the groundwater flow, salt transport and groundwater age distributions in coastal aquifer-aquitard systems examined in this study.This included both laboratory-and fieldscale simulations, whereby the model served to explore whether observations made at the laboratory scale were valid also for field-scale conditions.The relevant governing equations describing the flow and multiple-species solute transport in SUTRA-MS are given in the Supporting Information (i.e., Test S1).SUTRA-MS has been previously validated against laboratory experiments involving density-dependent flow and salt transport in coastal groundwater systems (Nguyen et al., 2020;Shen et al., 2016).SUTRA-MS uses Richards equation, with the unsaturated hydraulic conductivity and moisture content calculated using the formulae of van Genuchten (1980).A linear relationship between saltwater density (ρ s ) and salt concentration (C s ) was used, with ∂ρ/∂C set to a constant of 0.7143 (Lu et al., 2013).
Groundwater ages were simulated by introducing a second solute species (i.e., a time tracer) with a zero-order production rate of 0.0167 s 1 (i.e., 1 min 1 ), representing that the time tracer with a concentration of one is produced per unit minute.This allowed the time tracer to accumulate within the system, whereby the concentration (i.e., minute) represented the groundwater age.The time tracer was subject to the same advective and dispersive controls as the saltwater (Goode, 1996;Wilson & Gardner, 2006), except the tracer did not contribute to the water density (i.e., ∂ρ/∂C of the time tracer was set to zero).Any inflow, including seawater and freshwater, had an initial age of zero, while outflow occurred at the ambient age.This is the classic 'flow-in-flow-out' solute boundary condition, described by Voss and Souza (1987), that mimics the discharge of groundwater at the ambient solute concentration, whereas inflow occurs at a specified concentration, where inflow occurs at a specified concentration (e.g., zero).
The same types of boundary conditions were applied to the laboratory-scale and field-scale numerical models (Figures 2b and 2c).As surface processes such as recharge, evapotranspiration and surface runoff were neglected in both, the exposed surface (DE; see Figure 2) was treated as a no-flow boundary.The aquifer base (AF) was also a no-flow boundary.The landward boundary (E′F) was assigned a constant flux boundary, while the unsaturated part of the landward boundary (EE') was set to a no-flow condition.The aquifer-ocean interface (ABCD; i.e., red lines in Figure 2) was assigned a specified-pressure boundary condition based on the seawater density and the depth of inundation, with the depth of inundation depending on the sinusoidal tide, given by: where T is the tidal period [T], and t is time.A seepage face is expected to form during falling tides, resulting from the water table in aquifers decoupling from the sea level.This requires estimation of the exit point (i.e., the water table contact point with the shoreline; see Figure 1), which was obtained using the approach of Wilson and Gardner (2006) and Xin et al. (2009).They adopted the following strategy in simulating the seepage face during falling tide conditions: (a) boundary nodes above the sea level were set to atmospheric pressure if the node was saturated in the previous time step, (b) boundary nodes above the sea level were set to a no-flow condition if the node was unsaturated in the last time step, and (c) boundary nodes below sea level were assigned a specified-pressure condition reflecting the depth of inundation and the seawater density.The ocean boundary conditions for solutes involved a constant saltwater concentration assigned to nodes with inflow (into the aquifer), while a zero-concentration gradient was applied to nodes with outflow (into the sea).

Model Parameters and Simulations
All three of the laboratory experiments (Cases H, L1 and L2) and the field-scale scenarios (Case F) were simulated with SUTRA-MS.For the laboratory-scale model, the longitudinal dispersivity (α L ), transverse dispersivity (α T ), residual water saturation (S wres ), and van Genuchten (1980) parameters (including n and a) were adjusted via manual calibration by comparing the simulated salinity distributions to the plume images of the laboratory experiments.The resulting parameter values are given in Table 1.Low values of α L and α T were obtained, consistent with the narrow saltwater-freshwater mixing zone observed in the laboratory experiments (as shown in the experimental results; Section 3.1).The values of S wres , n and a for the coarse sand (i.e., the unsaturated zone only occurs in the top unconfined aquifer) were obtained through a comparison of the simulated and experimental extents of the USP.Tidal parameters and the inland constant flux used in the laboratory-scale model were based on experimental conditions.To guide the laboratory-and field-scale models, we introduced a dimensionless scale λ, which is defined as the ratio of the aquifer depth to the tidal propagation distance (Fang, Zheng, Wang, et al., 2022;Mo et al., 2021;Nielsen, 1990).
where ω is the tidal angular frequency (T 1 ), θ is the porosity ( ), and H is the aquifer thickness.

Water Resources Research
10.1029/2023WR035200 ZHANG ET AL.
The field-scale cross-sectional model was 400 m long and 48 m deep (below the top elevation of the onshore aquifer).The base case (i.e., Case F) adopted the parameters (K s , S wres , a and n) for sand suggested by Carsel and Parrish (1988) to represent the aquifer sediments (i.e., both the unconfined and semi-confined aquifers).The value for K s of the aquitard was set to 8.25 × 10 8 m s 1 (i.e., three orders of magnitude lower than the aquifer K s in Case F).The aquifer-to-aquitard hydraulic conductivity ratio (r′ K ) of 1,000 is consistent with the study of Hugman et al. (2015), who observed this value for a coastal aquifer-aquitard system in south Portugal.Knight et al. (2018) reviewed the value of r′ K for coastal aquifers from data given in prior literature, finding values ranging from about 100 to 10,000.
In the field-scale scenarios, the porosities of the aquifer and the aquitard were 0.45 and 0.48, respectively, following the numerical study of Zhang, Lu, et al. (2021).Dispersion parameters of α L , α T and D d (see Table 1) were set to the same values in all three layers (aquifer-aquitard-aquifer) consistent with the approach of Zhang, Lu, et al. (2021).Field-scale tides were defined by A = 1.5 m and T = 12 hr under mean sea level (H MSL ) at 45 m.A constant freshwater flux (Q) of 2.1 m 2 d 1 (volumetric flux per unit width of the aquifer) was applied to the landward boundary.This value was estimated in a field investigation of salinity distributions within a sand island aquifer (Moreton Island, Australia) of thickness ∼30 m (Robinson et al., 2006), and was adopted in subsequent modeling studies (e.g., Robinson, Li, & Barry, 2007;Xin et al., 2010).In layered-aquifer cases, inflows through the inland boundary to the upper and lower aquifers and the aquitard were set based on the transmissivity (K s multiplied by the initial saturated thickness) of the respective units.For example, the fractions of Q applied to the unconfined aquifer, aquitard and semi-confined aquifer of the field-scale model in Case F were approximately 4.25 × 10 1 , 6.45 × 10 5 and 5.75 × 10 1 , respectively.Note that the freshwater input applied to the aquitard along the inland constant-flux boundary is almost negligible.The calculation of λ for all laboratory-and field-scale models using Equation 2 is within the range for natural coastal aquifers (Robinson, Li, & Prommer, 2007).
A sensitivity analysis was undertaken using the field-scale model to examine the effects of tidal amplitude and the aquitard K s on the groundwater flow and salt transport processes in the aquifer-aquitard system.Other parameters were retained from the base case (Case F; see Table 1).This led to eight additional cases, as: (1) aquitard Discretization of the laboratory-scale and field-scale models represented a balance between model accuracy and run times.The laboratory-scale model domain was discretized into 103,329 nodes and 102,400 elements, while 116,725 nodes and 115,840 elements were used in the field-scale simulations.Voss and Souza (1987) recommended that numerical instabilities can be avoided in SUTRA if the model discretization meets the grid Péclet number criterion of P e ≈ ΔL/α L ≤ 4, with ΔL being the transport distance [L] between two sides of an element along the groundwater flow direction.Maximum P e values for the laboratory-scale and field-scale models were 2.5 and 1.2, respectively, both of which meet the Péclet stability criterion.

Simulation Diagnostics
Several diagnostic variables were adopted to characterize the numerical results of SUTRA-MS.Rates of leakage (i.e., the discharge per unit area), q [L T 1 ], across the aquitard were calculated based on the method of Post et al. (2007) taking into account variable-density effects: where h f is the equivalent freshwater head [L], z is the elevation [L], and ρ a denotes the average water density between the top and the bottom of an aquitard [M L 3 ].Positive values of q indicate upward leakage.A comparison between q values obtained using Equation 3 and fluxes reported by SUTRA-MS showed an excellent match (r 2 ∼ 0.99), as illustrated in Figure S1 in Supporting Information S1, suggesting that Equation 3 is effective in quantifying the leakage rate across the aquitard in the variable-density flow system.
The effects of tides were characterized using two metrics: (a) tide-induced reduction in the upward leakage (through the aquitard) (ΔQ′ L ), and (b) tide-induced change in the toe position (Δx′ t ).These are given as: where Q LS [L 2 T 1 ] and Q LT [L 2 T 1 ] are the net upward leakage through the aquitard for static (i.e., non-tidal) and tidal (i.e., the time-averaged leakage is used) conditions, respectively; x tS [L] and x tT [L] are the toe positions of the saltwater wedge for static and tidal (i.e., quasi-steady values) conditions, respectively.Figures 3a and 3d show that the homogeneous aquifer contained a classical, density-driven saltwater wedge, characterized by a slanting and narrow mixing zone.The presence of an aquitard in Cases L1-S and L2-S reduced the landward extent of saltwater relative to the homogeneous case.Two separate saltwater wedges were created in the unconfined and semi-confined aquifers, although the higher K s of the aquitard in Case L1-S (Figures 3b and  3e) led to a substantially smaller deviation from the saltwater distribution of the homogeneous case (Case H-S; Figures 3a and 3d).The aquitard also caused widening of the mixing zone, especially within and adjacent to the aquitard (e.g., Case L1-S; Figures 3b and 3e).Lu et al. (2013) also observed wider mixing zones in the lowpermeability layer that overlaid a high-permeability layer in their non-tidal analysis of seawater intrusion in layered coastal aquifers.

Salinity Distributions
Freshwater occurred within the semi-confined aquifer beneath saltwater in the overlying unconfined aquifer in Cases L1-S (Figures 3b and 3e) and L2-S (Figures 3c and 3f).This is the consequence of the restriction to upward leakage of fresh groundwater from the semi-confined aquifer, leading to a small amount of OFG in Case L1-S, while in Case L2-S, an extensive body of OFG occurred that reached the seaward boundary of the sand tank.Thus, the extent of OFG depended on the aquitard K s , consistent with previous studies (e.g., Bakker et al., 2017;Lu et al., 2013).The saltwater wedges appear to contain relatively uniform salinities in the semi-confined aquifer of the layeredaquifer cases (Figures 3b, 3c, 3e and 3f) and in the homogeneous case (Figures 3a and 3b).In contrast, the unconfined aquifers of the layered-aquifer cases had saltwater wedges containing irregular freshwater-saltwater mixing.This was apparent in the sand tank experiments as lighter-colored patches caused by vertical fingers of freshwater emanating from the underlying aquitard (Figures 3b and 3c).These classical features of unstable, mixed-convective flow are more clearly illustrated in the numerical modeling results (Figures 3e and 3f).Instabilities were apparent even in Case L2-S, where freshwater occurred beneath saltwater over only a small distance (Figures 3c and 3f).As illustrated in Figure S2 (see in Supporting Information S1), the distribution of the normalized pore-water velocity of Case L2-S demonstrates that the vertical freshwater flow crossed the aquitard into the overlying saltwater wedge, leading to the formation of unstable freshwater fingers.The freshwater fingers within the saltwater wedge of the unconfined aquifer (e.g., Figures 3e and 3f) were in continuous motion, while other elements of the experiment reached a quasi-steady-state condition.For example, the extent of the wedge, indicated by the toe location, was stable after 24 hr, even though the salinity distribution within the wedge continued to change.The occurrence of mixed-convective flow and instabilities within the upper aquifer is similar to the observations of Solórzano-Rivas et al. (2021), who used numerical simulation to evaluate simpler, vertical flow scenarios of a subsea, unconfined aquifer exposed to underlying freshwater.There is also some evidence of unstable mixed-convective flow in the simulation results of Lu et al. (2013) involving the lowest aquitard permeability, although they did not discuss this phenomenon.
Figure 4 compares the quasi-steady salinity distributions from laboratory experiments and corresponding numerical simulations for Cases H-T, L1-T and L2-T (tidal conditions).The geometry of the saltwater-freshwater interface (i.e., including the saltwater wedge and the USP) within sand-tank experiments is closely approximated by the contour from the model representing 50% seawater concentration.Thus, the strong model-experiment match obtained for the non-tidal conditions was also achieved with tides included.
Figure 4 shows that the tide-induced USP caused the freshwater to discharge to the sea deeper within the aquifer, at approximately the low-tide level in all three cases.This forced the saltwater wedges within the unconfined aquifers to reduce in size (relative to non-tidal conditions; compare Figures 3 and 4).The same behavior was observed by Kuan et al. (2012) in their investigation into the effects of tides on seawater intrusion in unconfined coastal aquifers.
The saltwater wedge within the semi-confined aquifers of Cases L1-T and L2-T also shifted seaward with the introduction of tides (compare Figures 4b, 4c, 3b and 3c), even though semi-confined aquifers did not show the presence of a USP (that caused saltwater wedge retreat in the unconfined aquifers).For example, the toe position of the saltwater wedge in the semi-confined aquifer retreated seaward by approximately 30% (from 1.68 m in Case L2-S (Figures 3c) to 1.18 m in Case L2-T (Figure 4c)).This is the first observation of a tide-induced reduction in the extent of saltwater in a semi-confined coastal aquifer.The seaward shift in the saltwater wedge within the semiconfined aquifer, with the addition of tides, led to an expansion of the OFG (compare Figures 4b, 4c, 3b and 3c).A more detailed analysis of the processes leading to OFG expansion with the addition of tides is included in Section 3.2.
A comparison between Figures 3 and 4 indicates that tidal fluctuations widened the mixing zones in both homogeneous and layered aquifer cases.Unconfined aquifers of the layered-aquifer cases retained the irregular freshwater-saltwater mixing within the interior of saltwater wedges with the addition of tides.That is, the enhanced freshwater-seawater mixing that accompanied tidal fluctuations did not eliminate the instabilities arising from the mixed-convective flow.
Figure 5 shows the field-scale salinity distributions for Cases F-S and F-T.The saltwater wedges in Case F-S reached steady-state conditions after 10,000 days.The corresponding pressure and solute distributions were adopted as initial conditions for Case F-T.The salinity distributions in Case F-T reached quasi-steady conditions after a further 1,500 days of tidal fluctuations.Unstable freshwater fingers, driven by upward leakage, formed within the saltwater wedge of the unconfined aquifer in both cases.The freshwater fingers continuously changed, whereby new fingers forming at the base of the wedge grew upwards and coalesced with existing, older freshwater fingers.The mixed-convective flow led to variations in the salinity of groundwater at the seafloor boundary.
The introduction of tides changed the toe position of the saltwater wedge, which shifted seaward in the semiconfined aquifer by approximately 30% (from 77.75 m in Case F-S to 54.45 m in Case F-T) due to the tide.At the time-averaged position of the toe, the salinity varied between 16 and 18 kg m 3 over the same frequency (determined from monitoring five successive tidal cycles).

Inter-Aquifer Leakage Rates
Figure 6 shows the spatial distributions of leakage rates across the aquitard obtained from the simulations of Cases F-S and F-T.Leakage rates from tidal simulations represent the sum over a tidal cycle (with positive values reflecting upward leakage) during which the direction of leakage fluctuates.Values were obtained after 3,000 tidal cycles (1,500 days) had occurred to obtain quasi-steady values.The total leakage through the aquitard in the non-tidal case (Case F-S) was 0.66 m 2 d 1 , while in the tidal case (Case F-T), the leakage (time-averaged) was 0.54 m 2 d 1 .Tides caused the leakage in the onshore region (i.e., the light orange shaded area; Figure 6) to decrease by 33%, while in the offshore part (i.e., the light green shaded area), leakage increased by 7% with the addition of tides.Surprisingly, tides caused only a minor change (1.9% increase) in the leakage through the sloping beach (i.e., the light pink shaded area), which contains the intertidal zone.
The results described above indicate that there was more leakage across the aquitard (by approximately 19%) in the non-tidal case.This is despite a fixed-flux condition at the inland boundary.The differences in the net upward leakage, caused by the addition of tides, led to changes in the rates of freshwater discharge from the semi-confined aquifer to the left, vertical boundary.That is, freshwater outflow to the offshore vertical boundary increased by approximately 20% from 0.57 m 2 d 1 in Case F-S to 0.683 m 2 d 1 in Case F-T.The increase of 0.113 m 2 d 1 in the freshwater discharge to the offshore boundary was equal to the reduction in upward leakage that occurred in the sloping beach and onshore regions due to the same fixed flux at the inland boundary in both the tidal and non-tidal Water Resources Research 10.1029/2023WR035200 cases.The enhanced flow toward the offshore boundary was accompanied by a shift in the time-averaged position of the wedge toe in the semi-confined aquifer, which retreated seaward by approximately 30% with the addition of tides (i.e., compare Figures 5a and 5d).
Figure 7 shows the discharge across the aquifer-ocean interface in the semi-confined and unconfined aquifers of Cases F-S and F-T.Rates of groundwater discharge to the ocean in tidal cases are integrated over the tidal cycle.Figure 7 differentiates the average fluxes occurring during periods of influx and efflux.Following the approach of Robinson, Li, and Barry (2007), the spatial patterns of influx and efflux across the sloping beach surface (i.e., in the unconfined aquifer) over the tidal cycle were divided into three distinct zones: (a) density-driven recirculation zone (DRZ), in which seawater flows into the aquifer, and high-salinity groundwater that includes recirculated and diluted seawater discharges to the sea, due in part to the seawater-freshwater density difference; (b) freshwater discharge zone (FDZ), where groundwater discharges to the sea with a salinity of less than 5% of seawater; (c) tidally driven recirculation zone (TRZ), which is defined by the time-averaged position of the USP.All three zones occurred within the unconfined aquifer, whereas only the FDZ and DRZ were present in the semi-confined aquifer.The outflow via the TRZ in the unconfined aquifer, caused by tides, created a substantial increase (63%) in the total water efflux (integrated across the sloping sea boundary) in Case F-T (Figure 7b) relative to Case F-S (Figure 7a).Tides also caused the total water influx and efflux within the DRZ of the unconfined aquifer to increase by 66% and 53%, respectively (i.e., from 0.41 m 2 d 1 to 0.53 m 2 d 1 in Case F-S to 0.68 m 2 d 1 and 0.81 m 2 d 1 in Case F-T, respectively).Water efflux exceeded the influx within the DRZ as a consequence of the seaward fluxes of fresh groundwater emanating from the inland boundary.Upward leakage through the aquitard into the base of the seawater wedge within the unconfined aquifer occurred at rates of 0.13 m 2 d 1 and 0.14 m 2 d 1 in Cases F-S and F-T, respectively.This roughly balanced the differences in water flux between the inflow and outflow through the DRZ of the sloping boundary in Cases F-S and F-T.It is expected that the formation of mixed influx-efflux patterns across the sloping sea boundary (Figures 7a and 7b) is attributable to upward leakage through the unconfined aquifer wedge.A comparison of Figures 7c and 7d shows the effects of tides on discharge across the vertical offshore boundary, whereby tides create time-varying influx and efflux within the semi-confined aquifer (Figure 7d), while non-tidal conditions produce substantially smaller influxes that are indetectable at the scale of the vertical axis in Figure 7c.For example, the total saltwater influx in the DRZ of the semi-confined aquifer increased by approximately 10 times from 0.067 m 2 d 1 in Case F-S to 0.763 m 2 d 1 in Case F-T.

Distributions of Hydraulic Head
Figure 8 shows the distributions in the hydraulic head (converted to equivalent freshwater head) obtained from the simulations of Cases F-S and F-T.The numerical results from Case F-T are presented as time-averaged hydraulic heads (Figure 8b), for high-tide conditions (Figure 8c), and for low-tide conditions (Figure 8d).As shown in Figures 8a and 8b, the time-averaged hydraulic heads decreased toward the sea in both unconfined and semiconfined aquifers of Cases F-S and F-T due to the seaward discharge of freshwater.The net upward leakage calculated by Equation 3(i.e., Figure 6) indicated upward hydraulic gradients (time-averaged) across the aquitard in both Cases F-S and F-T, which resulted from higher hydraulic heads (i.e., considering a buoyancy correction) in the semi-confined aquifer than in the unconfined aquifer.The introduction of tides led to higher time-averaged heads in both aquifers.The difference between non-tidal and time-averaged tidal heads is referred to as watertable overheight, which is the super-elevation of heads near the coast arising from tidal effects (Nielsen, 1990).Tide-induced increases in heads were greater in the unconfined aquifer than in the semi-confined aquifer, as assessed by a comparison of hydraulic heads between Cases F-S and F-T.For example, tidal fluctuations induced a time-averaged watertable overheight of 0.42 and 0.06 m at x = 170 m (determined from modelgenerated hydrographs at the base of each aquifer) in the unconfined and semi-confined aquifers, respectively.The higher watertable overheight in the unconfined aquifer caused a weaker upward hydraulic gradient and, consequently, a reduction in the net upward leakage from the underlying semi-confined aquifer.
Tidal propagation in the unconfined aquifer was more attenuated than in the semi-confined aquifer, as a consequence of the higher storativity of the unconfined aquifer (Ferris, 1952).The range of the groundwater tide (i.e., high tide minus low tide) was <0.01 m where x > 218 m, whereas in the semi-confined aquifer, the tidal range was around 0.28 m at the inland boundary (x = 400 m).Similar observations of differences in the tidal propagation within aquifers of a layered system were made by Chen and Jiao (1999).Differences between the groundwater fluctuation amplitudes of the unconfined and semi-confined aquifers led to tide-induced fluctuations in the direction of the hydraulic gradient (and therefore flow) between the two aquifers.

Particle Tracing and Groundwater Age
Figure 9 shows the tracer paths and particle residence times in Cases F-S and F-T.The time-averaged pore-water velocities were used to determine the movement trajectories and residence times of particles under tidal conditions, following the approach of Robinson, Li, and Barry (2007).The release positions of particles in Case F-S were designated as follows: A and B originated at the offshore vertical boundary; C, D, E and F originated at the sloping seaward boundary; G, H, I and J originated in the freshwater zone (i.e., landward of the extent of seawater in the aquifers).For Case F-T, an extra particle was released at K to capture the time-averaged seawater circulation paths in the USP.
Figure 9 shows that the characteristics of density-driven circulation in both cases varied significantly between the overlying and underlying saltwater wedges as a consequence of freshwater leakage into the unconfined aquifer.Streamlines A and B depict counter-clockwise circulation in the wedge of the semi-confined aquifer, while the wedge in the unconfined aquifer hosted multiple circulation cells containing counter-clockwise (i.e., streamlines C, E and F; Case F-S) and clockwise (i.e., streamline D; Case F-S) circulation.This led to shorter particle residence times for saltwater recirculation in the unconfined aquifer compared to those of the semi-confined aquifer (i.e., Figure 9a).The residence times of freshwater were larger in the semi-confined aquifer in both Cases F-S and F-T, by a factor of about 4-5 relative to the unconfined aquifer.The pathways of streamlines I and J demonstrate the two main pathways for deep, old freshwater to exit the semi-confined aquifer, namely through horizontal flow toward the left, sea boundary, and through upward leakage to the unconfined aquifer.
After tides were introduced (i.e., Figure 9b), the path lengths of streamlines A and B of 103 and 48 m (respectively) in Case F-S reduced to 46 and 17 m (respectively) in Case F-T, and consequently, the residence times of 25.0 and 12.2 years (respectively) for streamlines A and B in Case F-S decreased to 8 and 2.5 years (respectively) in Case F-T.That is, the mean flow velocities of streamlines A and B of 4.1 and 3.9 m yr 1 (respectively) in Case F-S increased to 5.7 and 6.8 m yr 1 (respectively) in Case F-T.Tides also caused the mean saltwater flow velocities in the unconfined aquifer wedge (i.e., streamlines C, D, E and F) to increase by approximately 29% relative to those obtained for Case F-S.Thus, tides increased rates of seawater circulation in the saltwater wedge of both the unconfined and semi-confined aquifers.The tide-induced seawater circulation (i.e., the USP), indicated by streamline K, caused the freshwater path lengths to deviate and lengthen (i.e., streamlines G and H increased in length from 99 to 105 m, respectively, in Case F-S, to 120 and 124 m, respectively, in Case F-T).This extended the freshwater residence times in Case F-T.Robinson, Li, and Barry (2007) made similar observations, showing that freshwater discharge bypasses the lower edge of the USP, leading to increased freshwater residence times in unconfined aquifers.Tidal fluctuations had the opposite (albeit weaker) effect on the freshwater residence times within the semi-confined aquifer, whereby streamlines I and J changed from 8.8 to 9.7 years (respectively) in Case F-S to 8.7 and 8.8 years (respectively) in Case F-T, while the path lengths increased from 160 to 250 m (respectively) in Case F-S to 180 and 252 m (respectively) in Case F-T.Thus, tides caused a modest increase in the mean flow velocities of freshwater in the semi-confined aquifer.
The calculations of residence times of particles, as given above, neglect the effects of mixing.Figure 10 shows distributions of groundwater ages obtained from the application of the advection-dispersion equation (Goode, 1996), thereby accounting for dispersive mixing effects, for Cases F-S and F-T.The oldest groundwater within the semi-confined aquifer occurred around the toe location of the saltwater wedges, as observed previously by Post, Vandenbohede, et al. (2013) for homogeneous coastal confined aquifers, without considering tides (Figure 10).The low-velocity seawater circulation flow contributed to the oldest groundwater age within the semi-confined aquifer.Unlike the semi-confined aquifer, however, Figure 10 shows that the older age of groundwater occurred within the interior of the saltwater wedge in the unconfined aquifers.By comparing Figures 5 and 10, it is evident that the extent of these zones with older groundwater age almost overlaps with the distributions of freshwater fingers within the saltwater wedge.This indicates upward leakage into the unconfined aquifer brought older freshwater into the overlying wedge, whereby the fresher groundwater within the unconfined aquifer wedge was older than the recirculating seawater.Note that the differences in groundwater age between fresher groundwater and recirculating seawater seem to be lower than their differences in residence times (i.e., Figure 9).This occurs because dispersive mixing effects allow the dilution of the older freshwater by the younger seawater.Older freshwater is also evident within the lower part of the unconfined aquifer, inland of the wedge, due to slow-moving leakage through the aquitard from the lower aquifer (Figure 10). Figure 10b shows that tidal effects significantly reduced groundwater ages within the mixing zone of the semi-confined aquifer, leading to a small difference in groundwater ages between the mixing zone and the overlying freshwater in the semi-confined aquifer.This result is consistent with the decrease in the residence times of recirculating seawater response to tides within the semi-confined aquifer, as observed in Figure 9.Moreover, tide-induced seawater recirculation led to young groundwater ages in the USP, which was encircled by older freshwater that was a mixture of upward leakage from the semi-confined aquifer and inflows from the right landward boundary.
Upward leakage into the unconfined aquifer brought older freshwater into the overlying wedge, whereby the fresher groundwater within the unconfined aquifer wedge was older than the recirculating seawater.Older freshwater is also evident within the lower part of the unconfined aquifer, inland of the wedge, due to slowmoving leakage through the aquitard from the lower aquifer.Tide-induced seawater recirculation led to young groundwater ages in the USP, which was encircled by older freshwater that was a mixture of upward leakage from the semi-confined aquifer and inflows from the right, landward boundary.

Aquitard Hydraulic Conductivity
Figures S3-S7 (see in Supporting Information S1) summarize model outputs from the sensitivity analysis of the aquifer-to-aquitard hydraulic conductivity ratio (r′ K ) in terms of the total leakage through the aquitard, the freshwater and saltwater discharge through the vertical sea boundary, the total SGD, and the toe position of wedges in the unconfined and semi-confined aquifers.Variations to r′ K to determine sensitivities were achieved through changes to the aquitard K s , whereas the aquifer K s remained unchanged for all cases.Thus, the sensitivity analysis targeted the aquitard K s .Figure S3 in Supporting Information S1 illustrates that the net upward leakage from the semi-confined aquifer declined monotonically with increasing log (r′ K ) in both non-tidal and tidal cases.This is the expected consequence of lowering the aquitard K s , causing the inter-aquifer flux to reduce.Figure S3 in Supporting Information S1 also shows that higher values of the tide-induced reduction in upward leakage (ΔQ′ L ) were obtained for lower values of the aquitard K s , indicating that tides have a relatively stronger influence on leakage rates in aquifers with lower aquitard K s .Figure S4 in Supporting Information S1 shows that increasing the value of log (r′ K ) greatly enhanced the freshwater discharge across the vertical sea boundary under both non-tidal and tidal conditions, but caused only a slight increase in the saltwater discharge across the vertical sea boundary (under both non-tidal and tidal conditions).This indicates that seawater recirculation rates were not particularly sensitive to the rate of freshwater discharge to the sea.As a consequence, the total SGD almost remained largely unchanged with the increase in log (r′ K ) from 2 to four (i.e., Figure S5 in Supporting Information S1).
The toe position of the wedge in the unconfined aquifer shifted slightly landward when the value of log (r′ K ) increased (i.e., Figure S6 in Supporting Information S1).This was caused by the reduction in the upward leakage and subsequently a drop in the horizontal flow in the unconfined aquifer that accompanied the lowering of the aquitard K s .Figure S7 in Supporting Information S1 shows that increasing log (r′ K ) shifted the toe position of the wedge seaward in the semi-confined aquifer (under both non-tidal and tidal conditions), due to the enhanced horizontal flow (and heads).The reduced extent of saltwater wedges within the addition of tides (Figures S6 and S7 in Supporting Information S1) was expected in both unconfined and semi-confined aquifers.The toe shift (Δx′ t ) in the unconfined aquifer varied monotonically with log (r′ K ) (i.e., Figure S6 in Supporting Information S1), while Δx′ t in the semi-confined aquifer first increased and then decreased when log (r′ K ) was increased from 2 to 4. This indicates that as semi-confined aquifers become more confined (i.e., smaller aquifer K s and higher log (r′ K )), tides have a reduced influence on the toe position in the semi-confined aquifer.

Tidal Amplitude
Figures S8-S12 in Supporting Information S1 summarize model outputs from the sensitivity analysis of the tidal amplitude (A) in terms of the total leakage through the aquitard, the freshwater and saltwater discharge through the vertical sea boundary, the total SGD, and the toe position of wedges in the unconfined and semi-confined aquifers.Figure S8 in Supporting Information S1 shows that increasing the value of A reduced the net upward leakage from the semi-confined aquifer.This is attributable to the associated increase in the tidal watertable overheight (with larger A) that is expected to occur mainly in the unconfined aquifer.The value of ΔQ′ L monotonically increased with increasing A (Figure S8 in Supporting Information S1); the expected consequence of stronger tidal forces on the upward leakage rate and the aforementioned watertable overheight effect.Figure S9 in Supporting Information S1 illustrates that increasing the value of A contributed to a significant increase in both freshwater and saltwater discharge through the vertical sea boundary, highlighting the important role of the tide in enhancing seawater recirculation, in addition to the reduced leakage through the aquitard for reasons given above.This was reflected also in higher values for the total SGD with a larger A, as illustrated in Figure S10 in Supporting Information S1.Figures S11 and S12 in Supporting Information S1 show that the toe position of saltwater wedges in both unconfined and semi-confined aquifers monotonically retreated seaward when A became larger.This is a direct consequence of the enhanced watertable overheight.

Discussion
Layered aquifers are ubiquitous in coastal areas, where inter-aquifer leakage could occur.This leakage flow would contribute to a large amount of water exchange and chemical flux between shallow and deep aquifer systems, affecting the distribution and quality of fresh groundwater resources.Understanding the mechanisms responsible for leakage-linked flow and transport is crucial for evaluating these processes.The findings from this study provide a preliminary understanding of inter-aquifer leakage effects on flow and salinity dynamics within coastal aquifer-aquitard systems.They have several implications for unveiling the complex characteristics of seawater intrusion, SGD and geochemical hotspots within layered coastal aquifers.

Implications for Seawater Intrusion
The occurrence of freshwater-seawater mixing zones is a typical characteristic of seawater intrusion in coastal aquifers.Many field investigations observed wide mixing zones, ranging from hundreds of meters to kilometers (Price et al., 2003;Xue et al., 1993).Although large transverse dispersivity, tidal fluctuations, heterogeneity in the hydraulic conductivity, and kinetic mass transfer effects were proposed to explain the formation of the wide mixing zone, the formation of a wide mixing zone in real-world coastal aquifers still remains the subject of debate (Lu et al., 2009;Robinson, Li, & Barry, 2007;Yu & Michael, 2022).This occurs because theoretical models incorporating these factors often fail to reproduce the wide mixing zones obtained from field measurements.The results from this study show that upward leakage from the semi-confined aquifer into the saltwater wedge created complex saltwater-freshwater mixing within the interiors of saltwater zones.These irregular mixing zones could be identified as parts of the broad mixing zone through field observations.Under such circumstances, a wide mixing zone can be obtained by interpolating the salinity of observation wells.This explanation provides an alternative plausible mechanism responsible for wide mixing zones observed in real-world coastal aquifers.
The seaward shift in the saltwater wedge of the semi-confined aquifer due to tides caused an expansion of OFG, which is poorly understood globally (Micallef et al., 2021;Post, Groen, et al., 2013;Yu & Michael, 2019).Tides also enhanced seawater recirculation, even in offshore aquifers considered in the current analysis.As many coastal aquifers are subjected to substantial pumping, further research is warranted into the interactions between pumping and tidal fluctuations within layered-aquifer systems, to build on the enhanced understanding of tides arising from the current study.For example, the pumping-induced shift in the extent of OFG has not been examined for semi-confined aquifers in which both pumping and tidal fluctuations are considered.Given that tides change the vertical leakage through aquitards and the OFG extent, it is likely that human activities that enhance vertical leakage, such as abandoned and poorly constructed wells, cause changes to the OFG in complex ways that depend on interactions between the tide, pumping and the aquitard-to-aquifer hydraulic conductivity ratio, amongst other factors.Faults may also disrupt the integrity of aquitards and create vertical flow conduits.The role of aquifer heterogeneity more generally on the interactions between tides and pumping-induced seawater intrusion in layered aquifers remains an area requiring additional research effort.This will require analyses that consider three-dimensional representations of coastal aquifers to capture the radial flow of wells and to incorporate alongshore flows, given the work of Knight et al. (2021) that shows links between onshore and offshore aquifers that include important alongshore pathways.

Implications for SGD
Previous studies of the impacts of various factors (e.g., tides, waves and varying inland input) on SGD are mostly based on single unconfined aquifers (Fang, Zheng, Wang, et al., 2022;Robinson et al., 2006;Xin et al., 2010).The changes in the water efflux caused by leakage and tides presented in this study reveal the distribution patterns of SGD through layered coastal aquifers.For example, the upward leakage into the overlying saltwater wedge creates complex density gradients that lead to the occurrence of both clockwise and counterclockwise seawater circulation cells (i.e., Figure 9).These normal and reverse seawater circulations form a mixed influx-efflux pattern across the aquifer-ocean interface.Moreover, the addition of tides results in higher rates of both fresh and saline SGD through semi-confined aquifers, in addition to the tidal-induced rise in the seawater recirculation rates through unconfined aquifers.These results also provide preliminary guidance to evaluate the SGD-borne chemical loads to the sea in real-world cases.For example, the mixed influx-efflux mode could create a scattered pattern of contaminant discharge.This is required to deploy the monitoring instruments at relatively low intervals to obtain contaminant efflux concentrations in field investigations.Moreover, the land-sourced chemical discharge through deep aquifers should be taken into account in the coastal nutrient budgets, though it is difficult to be quantified.

Implications for Biogeochemical Hotspots
The observations of mixed-convective processes within the wedges of the upper, unconfined aquifer in laboratory experiments and numerical simulations are likely to have implications for biogeochemical reactions within layered coastal aquifers.For example, greater mixing is expected between chemicals of terrestrial origin and the recirculating seawater.This may lead to enhanced aerobic biodegradation and denitrification within the saltwater wedge where the seawater is enriched in organic matter and oxygen.The shorter pathways of recirculating seawater that occur from the mixed-convective processes within the saltwater wedge and the resulting shorter residence times may create higher inputs of oxygen and dissolved organic matter, further enhancing aerobic biodegradation and denitrification (Anwar et al., 2014;Santos et al., 2008).Thus, the saltwater wedges of layered aquifers may serve as more reactive zones than those of homogeneous settings.The biogeochemical processes within the interior of the saltwater wedge, including those that host mixed-convective processes such as demonstrated here, is a topic worthy of further research effort (Heiss et al., 2017;Kim et al., 2020;Robinson et al., 2009).This research should consider the mixed saltwater influx-efflux patterns that are expected in the nearshore aquifer, given the dependency of seafloor ecosystems on nutrient loads (e.g., Cui et al., 2021;Lee et al., 2007).

Conclusions
The laboratory experiments and numerical modeling undertaken in this study lead to the following conclusions.
1. Unconfined aquifers within layered aquifer systems may comprise salinity distributions that are controlled by unstable, buoyancy-driven fingers of freshwater.The development of freshwater fingers created irregular mixing zones of terrestrial freshwater and seawater within the interior of the saltwater wedges considered in this study.This caused the occurrence of both clockwise and counterclockwise circulation within the overlying unconfined aquifer, leading to the formation of mixed saltwater influx-efflux patterns across the nearshore aquifer-ocean interface.2. The addition of tides reduced the net upward leakage from the semi-confined aquifer to the overlying unconfined aquifer.This increased the horizontal flow toward the sea, which in turn reduced the extent of the saltwater wedge in the semi-confined aquifer.Where the tidal amplitude was larger, this phenomenon had a stronger influence on the extent of seawater in the semi-confined aquifer.3. The higher rates of both fresh and saline submarine groundwater discharge, caused by tides, led to lower groundwater ages in the semi-confined aquifer.This effect was more pronounced when the aquitard hydraulic conductivity was lower or the tidal amplitude was higher.Lowering the aquitard hydraulic conductivity had only a weak influence on the seawater recirculation rate within the semi-confined aquifer, which was influenced more so by the addition of tides.
Notwithstanding these findings, the main limitation of the current study lies in the lack of support from field data.Moreover, several important knowledge gaps remain.For instance, most aquifer-aquitard systems worldwide are exposed to longer-term forcings such as spring-neap tides (Heiss & Michael, 2014;Robinson, Gibbes, et al., 2007) and seasonally varying inland groundwater inputs (Fang, Zheng, Guo, et al., 2022;Michael et al., 2005), resulting in complex, longer-term flow and salt transport characteristics.Although semi-confined aquifers are increasingly showing the impacts of groundwater pollution (Han & Currell, 2022;Kwon et al., 2021), the transport pathways and mixing dynamics of land-sourced contaminants in semi-confined aquifers remain poorly known.Most importantly, few studies have incorporated reactive solute transport to investigate the effects of inter-aquifer leakage on the geochemical environment of unconfined aquifers (Atteia et al., 2005;Love et al., 2000).Further research is required to address these gaps in knowledge and provide deeper insights into the flow and transport dynamics that arise in response to the complex factors that influence coastal aquifer-aquitard systems.

Figure 1 .
Figure 1.Conceptual diagram of a tidally influenced aquifer-aquitard system, consisting of an unconfined aquifer, an aquitard and a semi-confined aquifer.Inflow through the right-hand boundary occurs from inland aquifers.Leakage occurs across the aquitard according to head differences between the unconfined and semi-confined aquifers.The inset diagram shows the formation of a seepage face during falling-tide periods.

Figure 2 .
Figure 2. Schematic diagrams of: (a) Components of the laboratory setup; (b) Dimensions and boundary conditions adopted in numerical models of the laboratory-scale apparatus; (c) Dimensions and boundary conditions adopted in numerical models of a hypothetical field-scale, layered coastal aquifer.The bold red lines in (b) and (c) represent the seaward boundary conditions, while the bold blue lines identify the inland boundary conditions.Dimensions are given in meters as x and z coordinates in brackets.The orange-filled circles indicate the positions of vertices.MSL refers to mean sea level.

Figure 3
Figure3presents the steady-state salinity distributions from laboratory experiments and corresponding numerical simulations of Cases H-S, L1-S and L2-S (non-tidal conditions).The results show that the experiments produced offshore fresh groundwater (OFG) within the semi-confined aquifer.The contour representing 50% seawater concentration, extracted from the numerical modeling results, captures the extent of the saltwater wedge observed in the sand-tank experiments reasonably well in all cases.This is an indication of the reliability of the SUTRA-MS

Figure 3 .
Figure 3.Comparison of non-tidal experiments (left column): (a) Case H-S, (b) Case L1-S and (c) Case L2-S and the corresponding simulations (right column): (d) Case H-S, (e) Case L1-S and (f) Case L2-S.The black solid lines in (a), (b) and (c) represent the location of salinities that are 50% of the seawater concentration, while black dashed lines identify mean sea level.The white arrows show the simulated pore-water velocity field.The white lines in (e) and (f) indicate the extent of the aquitard.

Figure 4 .
Figure 4. Comparison of tidal experiments (left column): (a) Case H-T, (b) Case L1-T and (c) Case L2-T and the corresponding simulations (right column): (d) Case H-T, (e) Case L1-T and (f) Case L2-T.Simulation parameters are given in Table 1.The black solid lines in (a), (b) and (c) represent the location of salinities from the model that are 50% of the seawater concentration, while black dashed lines identify mean sea level.The white arrows show the simulated pore-water velocity field.The white lines in (e) and (f) indicate the extent of the aquitard.

Figure 5 .
Figure 5. Salinity distributions of Case F-S (left panel) and Case F-T (right panel) at 50-day intervals.The black dashed lines identify mean sea level.The red solid lines indicate the extent of the aquitard.The times of 10,000 and 11,500 days (where the latter represents 1,500 days of tidal fluctuations) for the figures in the top row represent steady-state and quasi-steady conditions for Cases F-S and F-T, respectively.

Figure 6 .
Figure 6.Spatial distributions of leakage rates (time-averaged values are given for the tidal case) across the aquitard obtained from the simulations of Cases F-S and F-T.The light green, pink and orange shaded areas represent the offshore part, sloping beach and onshore part of the system, respectively.

Figure 7 .
Figure 7.Comparison of water fluxes across the aquifer-ocean interfaces of the unconfined aquifer (sloping sea boundary; top row) and semi-confined aquifer (vertical sea boundary; bottom row) between Cases F-S (left panels) and F-T (right panels).Flow to the ocean through the horizontal seabed is given as the green-shaded region of Figure 6.The tidal discharge rates of Case F-T are time-averaged values.Here, DRZ, TRZ and FDZ represent the density-driven recirculation zone (i.e., saltwater wedge), the tidally driven recirculation zone (i.e., USP), and the freshwater discharge zone, respectively.

Figure 8 .
Figure 8. Spatial distributions of the hydraulic head (converted to equivalent freshwater head) for: (a) Case F-S; (b) Case F-T (time-averaged values); (c) Case F-T (hightide conditions); (d) Case F-T (low-tide conditions).The black solid lines represent 50% seawater concentration contours.The black dashed lines indicate mean sea level in (a) and (b), and high-tide and low-tide sea levels in (c) and (d), respectively.The red solid lines depict the extent of the aquitard.

Figure 9 .
Figure 9.Comparison of tracer paths and particle residence times (in years) between: (a) Case F-S and (b) Case F-T.The black dashed-dot lines represent 50% seawater concentration contours, while the black dashed lines indicate mean sea level.Capital black letters indicate the release positions of particles, while numbers in brackets are the corresponding particle residence times.The colors of the lines show the residence times along the tracer paths.The red solid lines depict the extent of the aquitard.

Figure 10 .
Figure 10.The distribution of groundwater ages for: (a) Case F-S and (b) Case F-T.The black solid and dashed lines represent 50% seawater concentration contours and mean sea level, respectively.The red solid lines depict the extent of the aquitard.

Table 1
Parameters Adopted in Laboratory-Scale and Field-Scale Models Note.Base case values are given except where numbers are in bold, which indicate values of K s and A that were adopted in sensitivity testing using the field-scale models.aCalibrated values based on the laboratory experiments.ZHANG ET AL.