Groundwater, salt and contaminant transport in a coastal aquifer system with a low‐permeability layer in the intertidal zone

Low‐permeability layer (LPL), formed by natural deposit or artificial reclamation and commonly found below the intertidal zone of coastal groundwater system, can retard the ingress of seawater and contaminants, and shorten the travel time of the land‐sourced contaminant to the marine environment compared with a homogenous sandy coastal aquifer. However, there is limited understanding on how an intertidal LPL, a condition occurred in a coastal aquifer at Moreton Bay, Australia, influences the groundwater and contaminant transport across the shallow beach aquifer system. We characterized the aquifer hydrological parameters, monitored the in situ groundwater heads, and constructed a 2‐D numerical model to analyses the cross‐shore hydrological processes in this stratified system. The calibrated model suggests that in the lower aquifer, the inland‐source fresh groundwater flowed horizontally towards the sea, upwelled along the freshwater–saltwater interface, and exited the aquifer at the shore below the LPL. Whereas in the upper aquifer, the tidally driven seawater circulation formed a barrier that prevented fresh groundwater from horizontal transport and discharge to the beach above the LPL, thereby directing its leakage to the lower aquifer. A contaminant represented by a conservative tracer was ‘released’ the upper aquifer in the model and results showed that the spreading extent of the contaminant plume, the maximum rate of contaminant discharge to the ocean, and its plume length decreased compared with a simulation case in a homogenous sandy aquifer. Sensitivity analysis was also conducted to investigate the characteristics of the LPL, including its continuity and hydraulic conductivity, which were found to vary along the beach at Moreton Bay. The result shows that with a lower hydraulic conductivity and continuous layer of LPL reduced the groundwater exchange and contaminant transport between upper and lower aquifer. The findings from the combined field and modelling investigations on the impact of an intertidal LPL on coastal aquifer systems highlight its significant implications to alter the groundwater and mass transport across the land–ocean interface.


| INTRODUCTION
Submarine groundwater discharge (SGD) is an important terrestrial nutrient and contaminant source for the nearshore marine environment (Burnett et al., 2003;Slomp & Van Cappellen, 2004;Uchiyama et al., 2000;Westbrook et al., 2005).SGD in shallow unconfined homogenous aquifers is well-studied, including the influence of various forcing conditions on fresh groundwater dynamics and salt transport, and their effect on land-sourced contaminant flux into the ocean (e.g., Robinson et al., 2018).Most studies of the SGD in unconfined aquifers considered the complex tidally driven, variable density groundwater hydrodynamics in shallow, nearshore sandy aquifers (e.g.Robinson et al., 2007b, Llopis-Albert and Pulido-Velazquez, 2015, Shen et al., 2019).These processes are usually of concern because of the impact of local coastal development on fresh groundwater resources and dependent ecosystems, for example, by sustainable pumping strategies (Baena-Ruiz et al., 2018;Renau-Pruñonosa et al., 2016), sea level rise (Ketabchi et al., 2016;Rasmussen et al., 2013), seawater intrusion (Baena-Ruiz et al., 2020;Colombani et al., 2016) as well as the pollution of coastal waters by land-derived contaminants (Robinson et al., 2009).
In single aquifer systems, fresh groundwater discharges to the sea over a dense seawater wedge (SW).An upper saline plume (USP) forms above the freshwater discharge zone (FDZ) due to tidal forcing on the beach face (e.g., as illustrated in Robinson et al., [2018]).Contaminant transport in these aquifers follows the groundwater through the FDZ (Robinson et al., 2007b).The development of the USP is controlled by beach morphology.For example, a flat beach slope can lead to an unstable USP, leading to the formation of salt fingers (Greskowiak, 2014).Shen et al. (2019) numerically examined the effect of an unstable USP on the transport of a land-derived contaminant in a homogenous aquifer, finding that the contaminant moved deeper further into the aquifer with an expansion of the contaminant plume compared with a stable USP case.
Recent studies show that a low-permeability layer (LPL) can often exist in shallow coastal aquifers several meters below the intertidal zone in low energy environments (Evans et al., 2020;Oki et al., 1998;Ratner-Narovlansky et al., 2020).Such LPL can be formed due to fluctuations in sea level over time and may widely exist in the aquifer system up to kilometres landwards (Bratton, 2010).In some cases, land reclamation in urban areas results in sand aquifers overlying pre-existing wetlands with lower permeability.In shallow sandy aquifers, the existence of an LPL could divide the aquifers into an upper unconfined aquifer and a lower semi-confined or confined aquifer.Field investigation of a shallow coastal aquifer (within 20 m) with the LPL located below the intertidal zone suggested that the hydraulic head in the upper unconfined sand aquifer is rarely affected by tidal oscillation, whereas the lower confined or semi-confined sand aquifer is strongly influenced by tidal oscillations (Carey et al., 2009).
Besides affecting the magnitude of the tidal oscillation, the presence of the LPL could modify the SGD to the ocean.For example, SGD from the lower aquifer can be enhanced due to leakage from the upper aquifer (Guo et al., 2012) as the groundwater flow through the LPL is predominantly in the vertical direction (Oki et al., 1998).
Numerical studies have been carried out to investigate the effect of LPL on SGD discharge.Evans and Wilson (2017) investigated a coastal aquifer with a LPL located in Georgia, USA and the results indicated the fresh SGD is one order of magnitude greater for the case with an LPL than for the case without an LPL.Riedel et al. (2010) numerically investigated the groundwater dynamics in a salt marsh with a silt layer near the low tide level, showing the groundwater from the upper sand aquifer mainly discharged to the surface water during the low tide and 90% of SGD occurred in the lower sand aquifer.However, none of the previous studies has analysed the effect of water density on groundwater flow, a factor was not discussed in these studies, which could significantly modify the groundwater flow in coastal aquifers (Robinson et al., 2018).Furthermore, the continuity of the LPL, which is usually uncertain in many field investigations, may also affect the groundwater dynamics in layered aquifer systems (Trglavcnik et al., 2018).
The presence of an LPL could also affect the salinity distribution and landward contaminant transport in the coastal aquifer system.For example, with the influence of the LPL, the toe of the SW in the lower confined aquifer moves seaward with the decrease of the LPL's permeability (Frind, 1982), whereas the freshwater-saltwater mixing zone widens (Lu et al., 2013;Shi et al., 2018).The toe location and the mixing zone in the lower confined aquifer are largely controlled by the contrast of hydraulic conductivities between the layers, layer thicknesses, inland freshwater head (Lu et al., 2013) and salt dispersivity (Solórzano-Rivas et al., 2019).Following the work of Shen et al. (2019), Zhang et al. (2021) numerically examined the effects of a LPL located below the intertidal zone on flow and solute transport, and the results show the LPL could retain the contaminant.The contaminant plume was restricted with the decrease of the hydraulic conductivity and positively related to the depth of the LPL (below the intertidal zone).
The aforementioned numerical studies all assumed an LPL elevation located below the intertidal zone.However, as provided by geological evidence, LPL often exists in intertidal zone of a low energy coastal area (Carey et al., 2009;Riedel et al., 2010) and can support high-diversity ecosystems (Murray et al., 2019;Paterson et al., 2019;Peterson, 1991).Therefore, it is important to understand how groundwater dynamics and contaminant transport in this layered aquifer system adjacent to a low energy intertidal zone.In our study of a bayside site (i.e., Moreton Bay, Australia), The LPL is located at an elevation within the intertidal zone (Figure 2), which may impact the exit point of land-sourced nutrients and contaminants, and could significantly affect the nearshore organisms.This study numerically explores the effect of this intertidal LPL on SGD and contaminant transport in a layered aquifer system, which is based on our site.We considered the effect of (1) hydraulic controls (tidal amplitude and landward head) and (2) LPL characteristics (hydraulic conductivity, continuity and outcrop location).Parameter ranges were based on field observations.

| Study site
The study site is a coastal area on the Western side of Moreton Bay, Queensland, Australia (Figure 1a). Figure 1b shows monitoring wells at the site with an inlet surrounded by intertidal sand flats located to the east.This zone supports high biodiversity (Lovelock et al., 2019).
At the study site, an upper aquifer (average thickness = 2.0 m) consists of reworked natural coarse sand from land reclamation.This sand layer is on top of a previous intertidal wetland which forms a mixed silty/sandy clay layer (average thickness = 0.8 m), regarded as an LPL in this study (Figure 2).Natural sand is observed below the LPL.Field data from the boreholes shows the LPL exists between GW08 and ML1, but its continuity and beach outcrop formation remain unknown.

| Field monitoring
Multi-level monitoring bores ML1A and ML1E with depths of 1.5 and 7.5 m below ground level, respectively, were used to monitor the groundwater head in the aquifer above and below the intertidal LPL near the shoreline (Figure 1b).The inland groundwater head was monitored from the borehole GW08, located 160 m landward of ML1, with its screen located at 3.5 m below ground level.Three (Figure 3a) were obtained from a nearby tidal station (Figure 1a).Tidal Principle lunar semi-diurnal (M 2 ) was the main constituent.The range of spring and neap tide was approximately 1.0 m and 2.2 m, respectively.

| Numerical model and boundary conditions
The hydrological processes in a 2-D cross-shore section of a coastal aquifer with an intertidal LPL (Figure 2), were simulated using SUTRA- MS, a variable-density groundwater flow and multiple-species solute transport model (Joseph & Sanford, 2004).The beach topography and nearshore bathymetry for the model domain were based on the site digital elevation model.The geometry, stratigraphy and boundary conditions of the conceptual model are illustrated in Figure 2. Groundwater flow at our site will not be strictly 2-dimensional, however, our simplified conceptualization is aimed at identifying the importance of the LPL location on SGD groundwater dynamics.The inland boundary, AH (Figure 2) was prescribed with a hydrostatic pressure.No flow was applied to the bottom boundary GH, the top boundary AB, and the seaward boundary FG.BCDEF are tidal boundaries.C, D and E represent three beach breaks.BC represents a steep seawall at the site with a slope of 0.5.CD, DE and EF represent the tidal sand flats with slopes of 0.02, 0.01 and 0, respectively.
During the falling tide, the groundwater table in the aquifer is higher than the sea level, forming a seepage face on the beach.The seepage was simulated by attaching a 'ghost node' with a relative atmospheric pressure of 0, to the seepage node where pore water pressure is positive, and the discharge through the seepage is calculated by the product of the pressure difference between the two nodes and a conductance of 0.01 (m s À1 ) (see Wilson & Gardner, 2006;Xin et al., 2009).The unsaturated boundary nodes above the tide and seepage face were assigned to be no-flow boundaries, by setting the conductance between the boundary nodes to the attached 'ghost node' to be zero.Seawater salinity was set to be 35 ppt.Despite occurring intermittently, rainfall on the modelled zone was ignored as the cumulative amount had a negligible impact on the nearshore hydrological processes during the monitored period, except the increasing inland head at the beginning of the period (Figure 3a).

| Model parameters and simulations
The hydraulic conductivity of the sandy aquifer (K sand ) was around 45-60 m day À1 using Hazen's equation (Hazen, 1892), which was estimated from the particle size distributions of the sand collected from the aquifer surface to a depth of 1.7 m at the site.A value of 50 m day À1 was adopted by comparing the simulated groundwater head with the observed data.The K LPL was set to be 0.05 m day À1 using the same method as the K sand .The LPL has a constant thickness of 0.8 m, with its surface located at 0.0 mAHD (Australian Height Datum).Given the LPL was recorded in all the bore logs, in the base simulation case, this thin clay layer was assumed to be continuous across the modelled domain between the monitoring boreholes.The aquifer properties applied in these cases are listed in Tables 1 and 2.
The simulation domain was discretized by 14,904 nodes and 15,184 elements.The mesh discretization is non-uniform with refinement around the aquifer-ocean interface where the USP and SW are present (Δx = 0.5 m, Δz = 0.1 m).The choice of dispersivity and discretization satisfies Péclet number criteria Pe ≈ ΔL/α L ≤4 (Joseph & Sanford, 2004).Grid independence tests were performed to ascertain the convergence and stability of the numerical solutions.The model was implemented with a time-step of 600 s.

| Numerical experiment
The model was initially run to a steady state with a constant inland head of 1.0 mAHD and without tide.The spring-neap tidal cycles were then introduced by the following equation and the model was continued to reach a quasi-steady state where salt distribution in the aquifer became stable.
where h T is the mean sea level (MSL) (mAHD), A i , ω i and φ i are the amplitude (m), frequency (rad h À1 ) and phase shift (rad), respectively, determined by fitting to the monitored tidal data (R 2 > 0.98).

| The effect of an intertidal LPL on submarine discharge
The base case is to further investigate the relationship between variation of boundary conditions and groundwater dynamics in the layered aquifers compared with previous numerical study (Zhang et al., 2021).The 91 days of observed data at the landward (GW08) and seaward boundaries were introduced to the model as boundary conditions.Tidal boundary for the base case was prescribed with a tidal signal using a time resolution of 10 min, monitored from a tidal station (27 19 00 S, 153 10 0 E) (Figure 1b).As for the LPL outcrop, the LPL extends at one beach slope break C (x = 180 m, Figure 2) located at the MSL in base case.This setting was identified by comparing the simulated groundwater head with field measurements (Table 2).
LI ET AL.

| The effect of an intertidal LPL on contaminant transport
A land-derived conservative tracer was injected in the aquifer based on the settings of the base to examine the potential effect of the LPL on contaminant transport in the layered aquifer system.The parameter settings were consistent with those in the base case, adopting a medium hydraulic conductivity (K LPL ) in the LPL (Case K-M).In addition, a relatively low K LPL of 5 Â 10 À4 m day À1 , and high K LPL of 50 m day À1 , which represents a homogeneous aquifer, were also implemented in the Case K-L and Case K-H, respectively.The representation of the outcrop for the three cases was set to be at x = 180 m same as the base case.Such assumption is later relaxed by implementing two breaks (Case K-D) in the LPL between the two boreholes (A and ML1) in the scenario for contaminant transport.
Two breaks were arranged at 60-70 m and 125-135 m in the x-direction.The settings of all the cases are also listed in Table 2.
The cases of contaminant transport lasted for 250 days.The landward boundary was set to be a hydrostatic pressure with water table located at 1.0 mAHD after the initial quasi-steady state.This was chosen not only because there is no long-term inland groundwater head data available, but also because implementing a constant inland head boundary allows the focus to be on assessing the effect of the LPL hydraulic conductivity on the contaminant transport.The tidal boundary was allocated by Equation ( 1).The landderived conservative tracer with an injection rate of 10 À3 m 3 day À1 , and free of salt, was released under the water table centered at x = 120 m and z = 0.2 m.The contaminant injection event lasted for 1 day with a total mass of 1000 kg released.The location of the injection point (Figure 2) was set above the LPL to ensure the contaminant plume would not be affected by the LPL in the initial conditions.3 | RESULTS

| Field monitoring results
The measured groundwater heads show the stratigraphic behaviour of the layered aquifers (Figure 3a).The LPL separated the sand aquifers into an upper unconfined aquifer and a lower confined aquifer.The water level in the landward bore (GW08) decreased from 1.36 to 0.61 mAHD during the three-month monitoring period.Despite such monotonic decrease, the groundwater level in the landward bore was always higher than that monitored in the intertidal bores.Rainfall occurred intermittently during the monitoring period, with a cumulative depth of about 144 mm.Although groundwater level seems to be responsive to rainfall, the rainfallinduced water level rise is too small to change the local hydrological processes.For example, an increase of about 35 mm was observed The LPL broke at x = 60-70 m and 125-135 m despite the form of its outcrop remaining the same as base case.
at bore GW08 after a rainfall of 23 mm day À1 occurred at day 42, which does not change the overall decrease of the groundwater table.Both the daily mean water level in upper (ML1A) and lower (ML1E) aquifers generally follow the same pattern as the daily mean sea level (DMSL).On the contrary, the gradual decrease of landward groundwater level had more pronounced impact to upper aquifer than the lower aquifer, particularly at the early stage of the monitoring period.
The hourly hydraulic heads in the upper (Figure 3b) and lower (Figure 3c) aquifers showed again tide affects the lower aquifer more significantly than the upper aquifer.This is consistent with the find-

| Groundwater dynamics
The simulated groundwater levels in the upper and lower aquifers using the parameters above for the base case generally agreed with the observed data (Figure 3). Figure 4 shows the simulated spatiotemporal salt distribution during the monitoring period.Despite the presence of an intertidal LPL, the development of salt fingers from the USP followed the same mechanism to that in homogeneous aquifers (Shen et al., 2019).Specifically, the FDZ was located between a SW formed at the freshwater-saltwater interface, and a USP is developed at the steep beach.Additionally, seawater infiltrated from the sand flat (i.e., mild beach), leading to the formation of salt fingers between the USP and SW.Separated periodically from the USP, the salt fingers moved seaward and later merged with the SW (see also Greskowiak, 2014;Shen et al., 2019).With the decrease of the inland hydraulic head, both the USP and SW expanded further landward, as a result, the FDZ becomes much narrower, and newly developed salt fingers become thinner after 60 days.Figure 7a shows the total efflux of SGD and influx of seawater across the beach during the monitoring period.Although the daily influx appears to be correlated with tidal oscillation, as higher tidal amplitude may cause more recharge to the aquifer (Riedel et al., 2010;Wilson & Morris, 2012), the daily influx pattern at this site seems not coincide with the tidal signals (Figure 7a).We found that the influx is correlated more strongly with DMSL, with a correlation of 0.52 (Figure 3c), than that with tide amplitude (R 2 < 0.1).This is likely because the influx, especially in the upper aquifer is only dependent on high tide level (R 2 < 0.49) rather than low tide level.The average influx from the upper aquifer was about 1.27 m 3 day À1 m À1 , accounting for a quarter of the influx (3.29 m 3 day À1 m À1 ) from the whole beach.The total efflux decreased gradually during the monitoring period because of the decline of the inland groundwater level (GW08) (Figure 3a), whereas the efflux from the upper aquifer was almost zero.The correlation (R 2 ) of the inland hydraulic head and the efflux was 0.76 (Figure 7b).However, no significant (R 2 < 0.1) correlation was found between the efflux and tidal amplitude, in contrast with the previous studies in coastal salt marshes (Riedel et al., 2010, Wilson & Morris, 2012).Furthermore, Wilson et al. (2015) suggested the efflux was negatively related to the DMSL, which was also not observed in this study.The difference may be because of coastal morphology that these early studies were conducted in salt marshes where the slope was relatively flat and the surface was inundated with seawater during high tide.However, our field site is above the high tide level with a sharp beach slope and thus the tidal effect was likely reduced.(x = 150 m), forming a preferential downward flow from the upper aquifer through the LPL.During the MHWS (Figure 8d), tidal force reversed the flow direction near the beach, a phenomenon that can be found when the inland fresh groundwater inflow is weak (Boufadel, 2000;Robinson et al., 2007b).Similarly, the reversed groundwater flow encountered the upcoming freshwater that flowed seaward, forming a freshwater upwelling pathway from the lower aquifer to the upper one, at the same location as the preferential downward flow channel.However, this phenomenon was not found during the MHWN (Figure 8b), suggesting there was little leakage from the lower aquifer when the tidal amplitude was small.

| The effect of the LPL on contaminant transport
The model was applied to investigate the transport of a conservative solute released in the upper aquifer.Multiple cases were implemented to assess the impacts of LPL hydraulic conductivity, and the presence of LPL breaks on the contaminant transport pathways, as shown in Figure 9.The salt distribution of Case K-M, where aquifer settings are the same as the base case, was discussed before (Figure 4).Comparing with Case K-M, (Figure 9b), Case K-H, which adopted the hydraulic conductivity of the sand aquifer to the LPL, produced an enlarged USP, a small FDZ and a retreated SW towards the sea (Figure 9a).On the contrary, with the decrease of K LPL , Case K-L produced an enlarged FDZ, a small USP and salt fingers (Figure 9c).
Despite that, the toe of the SW remained at the same location both in Case K-M and K-L, suggesting the hydraulic conductivity of the LPL affected little the density-driven circulation.
In Case K-M, the injected contaminant formed a plume that later from leaching into the lower aquifer in Case K-L (Figure 9c1).As a consequence, the plume could only spread horizontally and remained at a high concentration in the upper sandy aquifer.Part of the contaminant infiltrated and stagnated into the upper section of the LPL because of the low permeability, leaving a footprint along its transport trajectory (Figure 9c2).The interception of vertical flow boosted the horizontal transport of the contaminant, causing it arrived at the USP in a shorter time than the other two cases.With the block of the vertical pathway, the contaminant plumed was forced to discharge to the ocean by transporting directly across the center of the USP (Figure 9c3).As LPL may not be continuous in coastal aquifers (Evans & Wilson, 2017;Jankowski & Beck, 2000;Trglavcnik et al., 2018) Based on Case K-M, we further quantify the effect of K LPL on the contaminant transport in coastal aquifers (Case K-L and Case K-H).
We used normalized contaminant mass (defined as the fraction of the residual contaminant mass in the aquifer to its total mass initially released, M/M 0 ) (Figure 10a), temporal mass flux across the beach (Figure 10b), contaminant plume area (Figure 10c), the trajectory of the plume center (Figure 10d) and spatial total mass discharge across the beach (Figure 10e) for sensitivity analysis.During the first 73 days, the mass remained unchanged as predicted by the three cases, since the contaminant had not reached the beach interface (Figure 10a).
The contaminant started to discharge from the aquifer after 74, 75 and 90 days in Case K-H, K-L and K-M, respectively.The contaminant moved much slower in Case K-M than that in Case K-H, as the lower permeable LPL suppressed the impact of the tide-driven flow, which enhances contaminant discharge.Surprisingly, the contaminant moved faster in Case K-L than that of K-M, as the discharge pathway through the USP is much shorter than routing around it.The contaminant discharge slowed down after 170 days in Case K-L as large amount of mass was sealed in the LPL.By 250 days, less than 5% of the contaminant remained in the aquifer as estimated by Case K-M and K-H, much lower than 15% estimated by Case K-L.
The total contaminant discharge across the aquifer-ocean interface-exhibited periodic oscillating pattern (Figure 10b), which are influenced by the salt fingers developed along the tidal flat (Shen et al., 2019).For all the cases, the flux was absent for the first 50 days, increase to a peak during 100 and 150 days, and finally reduced to nearly zero.In Case K-H with a higher permeable LPL, the flux peak increased 55% and arrived 17 days earlier, compared with those found in Case K-M.Interestingly, in Case K-L with a lower permeable K LPL , the flux peak also increased and occurred earlier.Such phenomenon suggests that there is no monotonic correlation between the peak flux and the hydraulic conductivity of the LPL.After the flux peak in 117 days, Case K-M produced greater contaminant flux than the other two cases, and the outflux of the three cases reduced to almost the same rate after 250 days.To quantify the LPL effect on the contaminant plume size, the plume area over time was presented in The LPL also altered the plume pathway as evident by the trajectory of the plume centroid (Figure 10d).For Case K-H and K-M, the plume centroid moved passed the LPL following the same pathway.
Subsequently, the plume trajectory estimated by Case K-M redirected and moved seaward and upward, before the arrival at the beach.The  10e).For Case K-M, the peak of total flux was observed at a point between C and D, despite the fluxes at the breaks remaining relatively higher than those in the surrounding area.For Case K-L, the majority of contaminant exited the aquifer at C (MSL) consistent with the finding in Figure 10d.

| The effect of the form of the LPL outcrop
We explored the effect of the of LPL outcrop form (truncated, extended or beyond) on the groundwater transport with the model.
The results for the Case L-165, where the LPL does not intersect the beach (truncated at x = 165 m), indicates groundwater fluctuation in the upper aquifer (Figure 11a) was enhanced because of the outcrop break, whereas the peak of the water head in the lower aquifer at ML1E decreased approximately by 0.38 m, as compared with the base case (Figure 11b), the peak of the hydraulic head in the lower aquifer estimated by the Case L-165 appeared to be 1 h later than that estimated by the base case, suggesting tidal forcing was damped when LPL outcrop was absent (Li & Jiao, 2001) were alternately located in L-165 that is, influx followed by efflux and influx.This suggests the LPL outcrop to the beach surface plays an important role in controlling the distribution of SGD.
Despite the form of the LPL outcrop would influence the local exchange between surface water and groundwater, the net fluxes across the beach by the three cases, were almost the same (between 1.28 and 1.37 m 3 m À1 day À1 , Table 3).This consistent net fluxes produced by different outcrop settings confirmed the finding by (Evans & Wilson, (2017)) that the LPL has little effect on the net flux.et al., 2019), while allowed the USP to be located landward in the upper aquifer (Figure 13b).The impeded USP by the LPL outcrop made the FDZ pathway more evident (Case L-195, Figure 13b).

| DISCUSSION
Groundwater flow and contaminant transport processes in layered aquifer systems have been the subject of recent research (Evans et al., 2020;Pu et al., 2020;Tal et al., 2017;Trglavcnik et al., 2018).
However, previous studies focused on the coastal aquifer system with the LPL elevation located below the intertidal zone, and no study investigated land-derived contaminant transport in coastal aquifers with an intertidal LPL.The field monitoring and numerical results presented herein suggest complex hydrodynamics and contaminant transport in the intertidal layered aquifer system with several outcomes.
Previous research suggests that the SGD of Moreton Bay ranges from about 2.9 Â 10 À2 to 5.9 Â 10 À2 m 3 m À2 day À1 (Stewart et al., 2015).Our study found the average SGD is 4.5 Â 10 À2 m 3 m À2 day À1 ranging from 4.3 to 7.0 Â 10 À2 m 3 m À2 day À1 at the study site.In addition, consistent with the findings by Evans and Wilson (2017), this study confirmed the presence of the LPL has little effect on net discharge of groundwater to the ocean in the shoreline.However, without the retardation effect of LPL, the water and solute exchange at the interface could drastically (maximum 70%) increase, and thus the discharge of land- derived nutrient to the sea would be a significant source (Gibbes, Robinson, Carey, et al., 2008).Furthermore, groundwater dynamics in the confined aquifer are sensitive to the tidal signal propagation (Li et al., 2004;Trglavcnik et al., 2018) and mean sea level, highlighting the vulnerability of coastal groundwater resources to sea level rise as a result of global warming (Werner & Simmons, 2009).
Consistent with previous study (Zhang et al., 2021), the size of the USP and salt fingers decreased with K LPL decreasing (Figure 9).
With the development of salt fingers being restricted, the FDZ expanded along the beach, however, the FDZ was also restricted by the LPL outcrop (Case L-195, Figure 12).Zhang et al. (2021) indicated the contaminant rarely intruded the LPL and exited the aquifer faster with lower K LPL located below the intertidal zone compared with that with a higher K LPL .The exit point of the contaminant was below the MSL shoreline despite the change of the K LPL .However, when the LPL is located in the intertidal zone, the contaminant could exit the aquifer slower with the middle K LPL (K-M) compared with that with the high K LPL (K-H) and low K LPL (K-L).This finding has an implication that the contaminant distribution in the layered aquifer is strongly dependent on the hydraulic conductivity and location of the LPL.For example, when the LPL located in the intertidal zone is discontinuous or composed of sandy loam or silt with relatively high hydraulic conductivity, the contaminant released in the upper aquifer may spread downward to the lower aquifer.Such spreading allows the coastal aquifer to act as a contaminant buffering zone, as contaminant concentration in the discharging groundwater is reduced, and the contaminant discharge time to the marine system is increased.On the contrary, if the layer is continuous and mainly composed of low permeable materials, such as clay, the contaminant transport may only occur in the upper aquifer.Although the plume spread is reduced, the contaminant is more concentrated, with a reduced residence time.Results from this study further complemented and expanded the previous research by Zhang et al. (2021).The results concluded that the geological units with low permeability in the intertidal zone may significantly affect the contaminant spread in the aquifer.In addition, it is worth nothing that other parameters may affect the model results.For example, the toe of SW could move further seaward with the increase of contrast of hydraulic conductivities between each layer, inland freshwater head and salt dispersivity (Lu et al., 2013;Solórzano-Rivas et al., 2019).The LPL thickness may also play a key role in contaminant transport as the residence time and plume size could increase with the increase of LPL thickness (Zhang et al., 2021), especially if the LPL thickness fully or partly covered the tidal range.
Although Most studies, including this one, focusing on groundwater dynamics and contaminant transport layered aquifer systems, assumed each aquifer is homogeneous and isotropic (Evans & Wilson, 2017;Lu et al., 2013;Riedel et al., 2010;Trglavcnik et al., 2018).However, aquifers are usually heterogeneous in nature, which leads to complex SGD and contaminant transport (Geng et al., 2020;Kreyns et al., 2020).

| CONCLUDING REMARKS
Through field monitoring and numerical simulation, this study investigated the effects of an intertidal LPL on the groundwater dynamics, salt distribution, and land-derived contaminant transport in a layered coastal aquifer system.The extensive LPL divides a sandy aquifer into an upper unconfined aquifer and a lower confined aquifer with the whole system exposed to a tidal boundary (i.e., ocean).The lower aquifer becomes unconfined where it outcrops.Major findings include: • Field data show the tidal oscillations observed in the lower aquifer 20 m inland are consistent with the magnitude expected for a confined aquifer.However, the daily average hydraulic head in both the upper and lower aquifers reflected the DMSL.Comparison with the field hydraulic head data, the 2D modelling indicates the K LPL can be three orders of magnitude smaller than K sand .For this case, rather than directly discharging into the sea, the groundwater in the upper aquifer mainly flowed into the lower aquifer through the LPL before discharging through the two breaks in slope at C and D. With the decrease of K LPL , the development of USP and salt fingers was restricted with expanding of the FDZ.
• The form of the LPL outcrop affected the groundwater exchange between the upper and lower aquifers near the shoreline.Compared with the case of the LPL extending to, or beyond, the shoreline at the MSL elevation, simulations for a truncated LPL (Case L-165) showed that the influx and efflux was significantly larger across the tidal boundary.However, the net flux across the shoreline remained similar in the three cases.Furthermore, with the LPL outcrop beyond the shoreline at the MSL elevation, the development of the USP and salt fingers was limited along the interface.
• Simulations of a conservative contaminant transport from a source in the upper aquifer considered the effect on groundwater flow (and solute transport) of the magnitude of K LPL as well as LPL discontinuity.The spread of the contaminant plume and its pathway were reduced with a decrease of K LPL .For the case where K LPL divided by K sand equals 10 À3 (Case K-M), the contaminant migrated downward in the lower aquifer through the LPL, before approaching the USP.The residence time of the contaminant was prolonged with flux peak being less than that in the homogenous aquifer consistent with previous study with the LPL located below the intertidal zone (Zhang et al., 2021).However, for the where K LPL divided by K sand equals 10 À5 (Case K-L), the contaminant primarily moved in the upper aquifer with some effectively 'trapped' in the LPL.Consequently, the residence time decreased with the increase of flux peak.Thus, the LPL may not be able to reduce the contaminant risk to the ocean, once the hydraulic conductivity of the LPL is lower than a certain threshold.
• The LPL break between the solute source and shoreline could facilitates the solute to move into the lower aquifer once the contaminant passed the LPL break.
This study only considered the LPL located in the intertidal zone at a 'nearshore' scale.In addition, we neglected aquifer heterogeneity and elevation variation of the LPL, which may affect the interaction of fresh groundwater and seawater.Future studies should aim to reveal the effects of these factors on coastal groundwater dynamics and contaminant transport.
pressure loggers (Level Troll 500, In Situ Inc., USA) were installed in these boreholes to record the hydraulic heads every 5 min over a three-month period.Air pressure and rainfall data were obtained from a weather station, located about 3.3 km away from the monitoring site.Data were collected from March to June 2018.Tide data F I G U R E 2 Model geometry and boundary conditions.The beach slopes of BC, CD and DE are 0.5, 0.02 and 0.01, respectively.The location of GW08 was set to be the origin of the x-axis and the mean sea level (MSL) was set to be the origin of the y-axis.The grey frame shows the other two types of the low-permeability layer outcrop in the calibration stage.F I G U R E 1 (a) The location of the study site in Moreton Bay, QLD, Australia.The blue star shows the location of a tidal station.(b) Overview of the study site during a low-tide period when offshore tide flats are visible.Monitoring wells are shown by the black circles.ML1A and ML1E represent multi-level wells to monitor the pore water pressure in the upper and lower sandy aquifer.

F
I G U R E 3 (a) Observed daily rainfall, observed daily mean sea level (DMSL), and observed and simulated daily averaged groundwater heads at the monitoring bores locations.Observed and simulated hourly-averaged groundwater heads of (b) bore ML1A, (c) bore ML1E, and sea water level obtained from Maritime Safety Queensland during 50 to 65 days.Comparison of observed and simulated water heads of (d) ML1A (d) and (e) ML1E, where R 2 equals 0.54 and 0.98, respectively.(f) Elevation of LPL (grey rectangle) with sea water level obtained from Maritime Safety Queensland.Elapsed days are calculated from 27 March 2018.The orange dots, purple triangles, and green dashed-dotted lines represent observed water heads of ML1A, ML1B, and GW08, respectively.The blue dashed-dotted lines represent DMSL in (a), and tide in (c) and (f) The red lines and blue lines represent simulated water heads of ML1A and ML1B, respectively.
entry suction parameter, a (m À1 ) 14.5 a 2.5.3 | The effect of the LPL outcrop on hydrodynamicsTwo more cases aimed to investigate the shape of LPL outcrop by comparing with the results from base case, as it is unknown during the field investigation.The two cases adopted the same aquifer hydrological properties, but with different formations of LPL outcrop (truncated, and beyond at the MSL shoreline).According to the location of the LPL outcrop, the base case can also be named Case L-180.The LPL formation of the second case was truncated at x = 165 m, landward to that of base case.This case is therefore named Case L-165.The third case allows the LPL to be extended beyond the beach at x = 195 m, forming an outcrop below the MSL.Similarly, this case is named Case L-195.
ings ofCarey et al. (2009) andTrglavcnik et al. (2018), indicating a confined aquifer because of the LPL is present to reduce the tidal damping.During high tide, the hydraulic head of the lower aquifer was at the same height as the sea water level, and the head dropped following the falling tide.During the low tide, the head remained constant at the MSL, suggesting the drainage slowed in the lower aquifer(Abdollahi-Nasab et al., 2019).

F
I G U R E 4 Daily average spatiotemporal variations of salinity distributions (daily average) for the base case.The yellow dashed line represents the mean groundwater level during the tidal cycle.The red and blue dashed lines represent the high tide (HT) and low tide (LT) during the day, respectively.The grey lines represent the extent of the low-permeability layer.

Figure 5
Figure5shows the tidally averaged groundwater level, salt distribution and the pathways of particles released from either landward or seaward boundaries.Landward fresh groundwater transported horizontally towards the sea, upwelled following the SW, and exited the aquifer at the beach below the LPL.The landward freshwater in the upper aquifer flowed downward into the lower aquifer, and eventually discharged into the sea via FDZ.The presence of the LPL significantly lengthened the travel time of the freshwater in the upper aquifer (about 1705 days) as well as the period of the tide-induced circulation in the USP (about 1738 days).Furthermore, both anticlockwise and clockwise density-driven circulations were observed along the sand flat, possibly due to the beach break near the mean low water neaps (MLWN) and the flat slope of the offshore beach(Zhang et al., 2017).Despite the variation of tidal amplitude during the spring-neap tide cycle, the exit point of the FDZ remained

F
I G U R E 5 Simulated average particle pathways under a springneap tide cycle.The numbers show the travel time of particles in days.The travel time is calculated based on the length of the pathway and the velocities along the line.The yellow dashed line represents the mean groundwater level during the tidal cycle.The red dotted lines represent mean high-water springs (MHWS) and mean low-water springs (MLWS).The blue dashed-dotted lines represent mean highwater neaps (MHWN) and mean low-water neaps (MLWN), respectively.The grey lines represent the extent of the lowpermeability layer.F I G U R E 6 Efflux and influx across the beach during spring tide (spring), neap tide (neap) and one spring-neap tide cycle (average).Positive value represents efflux.The grey line represents the beach elevation with two beach slope breaks (x = 180 and 195 m).The grey dashed lines represent the location of the beach breaks.

Figure 8
Figure8shows the simulated flow velocities between the upper and lower aquifers near the USP during a spring-neap tide.As seawater receded during the MLWN, groundwater in the lower aquifer flowed seaward and slightly downward, and the exchange between the two aquifers appeared to be very weak (Figure8a).In contrast, during

F
I G U R E 9 Spatiotemporal transport of land-derived contaminant released in the upper aquifer for Cases K-L, K-M, K-H and K-D.Red solid line represents salinity contour line (10 ppt).Yellow dashed line represents the average water table under a spring-neap tidal cycle.The red and blue dashed lines represent the high tide (HT) and low tide (LT) during the day, respectively.The grey lines represent the location of the lowpermeability layer.
moved horizontally in the upper aquifer and expanded into the lower aquifer through the LPL (Figure 9b1-12b2).The plume in the upper aquifer became stagnant around the 150 m point at the x-axis (Figure 9b3), as it was intercepted by the USP, which is consistent with the groundwater flow patterns shown in Figure 8. Consequently, the plume moved into the lower aquifer before its discharge along the FDZ.As a comparison, the contaminant plume produced by Case K-H moved much faster both horizontally and vertically in the upper section of the homogeneous aquifer (Figure 9a1-a2).However, with the K LPL decreasing to a lower value, while still in the range of the clay soils measured at the study site, the LPL prevented the contaminant F I G U R E 1 0 Temporal variations in (a) the normalized total mass (M/M 0 ), (b) mass flux through the beach, (c) plume spreading area, and (d) trajectory of the plume centroid, (e) spatial distribution of total contaminant discharge across the beach in Case K-H (black dashed line), K-M (black solid line) and K-L (grey solid line).The red solid line represents the beach slope.The blue dot-dashed line represents the low-permeability layer.C and D in (d) and (e) represent two beach breaks located at x = 180 m and 195 m, respectively.
, we hereby investigated the contaminant transport in the aquifers with two LPL breaks at 60-70 m and 125-135 m at the cross-shore axis (Case K-D, Figure 9d).The plume by Case K-D escaped from the upper aquifer to the lower one through the breaking cavity and infiltrated to the same depth as that achieved by Case K-H with a homogenous aquifer after 75 days.Consequently, less contaminant remained in the upper aquifer downstream to the break.The remaining contaminant in the upper aquifer still moved into the lower aquifer through the LPL at the same location as the plume achieved in Case K-M (x = 150 m).

Figure 10c .
Figure 10c.Similar to the Mass flux across the beach, the plume area grew from the release point in the beginning, expanded to a greatest extend occurred around 100 days, and shrunk thereafter as the contaminant discharged from the aquifer.Increasing the hydraulic conductivity of LPL from the base case (Case K-M) resulted in a larger plume area by 50 m 2 at its peak, and faster contraction following the peak, as shown the Case K-H.Decreasing the hydraulic conductivity of LPL reduces the peak plume area, and locked in part of the contaminant in the LPL, leading to a greater plume size at the end of the simulation in Case K-L.
increase of hydraulic conductivity in LPL allows the plume to transport deeper in the lower aquifer before the redirection.Such later redirection in Case K-H caused the discharging point seaward to that of Case K-M.Different to the previous two cases, Case K-L predicted the trajectory remained only in the upper aquifer and the intersected with the beach at the slope break C.The spatial distribution of contaminant discharge along the aquifer-ocean interface showed notably different patterns between Case K-H, K-M and K-L.Without the LPL (Case K-H), contaminant exited mainly in the two beach slope breaks (C and D) with the total flux of 55.2 and 72.5 kg, respectively (Figure . The groundwater fluctuation in the upper aquifer estimated by Case L-195, where LPL outcrop presents at the beach, was similar as that by the base case.However, the average hydraulic head of the lower aquifer was 0.14 m lower than that by the base case, likely due to the reduced leakage through the LPL in Case L-195.The two peaks of efflux across the aquifer-ocean interface over a spring-neap cycle estimated by Case L-195 appeared near the two beach breaks, respectively (x = 180 and 190 m, Figure 12), consistent with the flux estimation by the base case (Figure 6).Although the peak effluxes at the breaks by Case L-195 were greater than those by the base case.No exchange of surface water and groundwater occurred at the beach covered by the LPL (between x = 180 m and x = 195 m) as predicted by Case L-195, suggesting the LPL outcrop restricted the development of seepage between the breaks.Both the peaks of efflux and influx at the landward beach break (x = 180) given by Case L-165 were greater than those by base case.Furthermore, the efflux and influx along the tidal flat (from 195 m in the x-axis) Despite the similar net fluxes, the presence of LPL outcrop on the beach (Case L-195) reduced both the efflux and influx almost by half, especially in the lower aquifer.The reduced water influx because of the existence of LPL outcrop (Case L-195) brought in less salts into the beach, resulting in less total salt mass in both USP and SW, and a much smaller size of USP (Figure 13b).Comparing the USP formed by Case L-165 and that by Case L-195, a continuous LPL with the outcrop beyond the beach at MSL (Case L-195) restricted the salt fingers to expand vertically from the USP to the deeper aquifer (Shen

F
I G U R E 1 1 Temporal simulated water heads in the upper and lower aquifer in the cases when the length of the low-permeability layer bottom (L LPL ) was 165 m (L-165), 180 m (Base), and 195 m (L-195), respectively.F I G U R E 1 2 Efflux and influx across the aquifer-ocean interface by Cases LPL-165 (blue lines) and LPL-195 (red lines).Positive value represents efflux and negative value represents influx.The grey line represents the beach elevation with two beach slope breaks at x = 180 m and 195 m.The grey dashed lines represent the location of the beach breaks.

F
I G U R E 1 3 Average salt distribution by Cases L-165 (a) and L-195 (b) under a spring-neap tidal cycle.The yellow dashed line is the average water table under a spring-neap tidal cycle.The grey lines represent the location of the low-permeability layer.
the aforementioned studies of SGD have investigated the transport process in a representative coastal aquifer with a 2-D crossshore vertical model, Zhang et al. (2016) compared salinity distributions at a beach aquifer for 2-D and 3-D models.The saltwater intruded further landward in the 3-D model than in the 2-D model, due to the alongshore flow and a nearby tidal creek.Groundwater flow at the site is likely to be affected by a more geometrically complex tidal boundary.Thus, to better understand the groundwater dynamics at the site scale, a 3-D model could be used to investigate this coastal area.
Parameters for sensitivity analysis.
T A B L E 2 a Comparison of the LPL continuity effect on the water fluxes across the beach, and salt mass in the aquifer (total), USP and SW.